% Optimalizace - Diplomová práce - Růžička Pavel

clc
clear all
close all

% Vstupni data
n=8197;
T0=14.715;

% Načtení dat z Excelového souboru do matice
filename = 'Data_optimalizace_vstup.xlsx';
sheet = 'List1'; 
vstupni_data = readmatrix(filename, 'Sheet', sheet);

alfa = vstupni_data(:, 1);
K = vstupni_data(:, 2);
thetatwist = vstupni_data(:, 3);
T = vstupni_data(:,4);
Q = vstupni_data(:,5);

alfa=alfa';
K=K';
thetatwist=thetatwist';
T=T';
Q=Q';

P = (2*pi*n/60)*Q;

%Matice proměnných
x1 = K;
x2 = alfa;
x3 = thetatwist;
X = [x1', x2', x3'];

% regresní model
regrese_T = fitlm(X, T, 'quadratic'); % kvadraticka regrese - tah
regrese_P = fitlm(X, P, 'quadratic'); % kvadraticka regrese - vykon

% Koeficienty regrese
disp('Koeficienty - tah:');
disp(regrese_T.Coefficients);
disp('Koeficienty - výkon:');
disp(regrese_P.Coefficients);

% Stanovené koeficienty z regresní analýzy
a0 = regrese_T.Coefficients.Estimate(1);
a_x1 = regrese_T.Coefficients.Estimate(2);
a_x2 = regrese_T.Coefficients.Estimate(3);
a_x3 = regrese_T.Coefficients.Estimate(4);
a_x1x2 = regrese_T.Coefficients.Estimate(5);
a_x1x3 = regrese_T.Coefficients.Estimate(6);
a_x2x3 = regrese_T.Coefficients.Estimate(7);
a_x1sq = regrese_T.Coefficients.Estimate(8);
a_x2sq = regrese_T.Coefficients.Estimate(9);
a_x3sq = regrese_T.Coefficients.Estimate(10); 

% Stanovené koeficienty z regresní analýzy
b0 = regrese_P.Coefficients.Estimate(1);
b_x1 = regrese_P.Coefficients.Estimate(2);
b_x2 = regrese_P.Coefficients.Estimate(3);
b_x3 = regrese_P.Coefficients.Estimate(4);
b_x1x2 = regrese_P.Coefficients.Estimate(5);
b_x1x3 = regrese_P.Coefficients.Estimate(6);
b_x2x3 = regrese_P.Coefficients.Estimate(7);
b_x1sq = regrese_P.Coefficients.Estimate(8); 
b_x2sq = regrese_P.Coefficients.Estimate(9);
b_x3sq = regrese_P.Coefficients.Estimate(10); 


% Generování rovnice
T_rov = sprintf('T = %.4f + %.4fx1 + %.4fx2 + %.4fx3 + %.4fx1x2 + %.4fx1x3 + %.4fx2x3  + %.4fx1^2 + %.4fx2^2 + %.4fx3^2', ...
                    a0, a_x1, a_x2, a_x3, a_x1x2, a_x1x3, a_x2x3, a_x1sq, a_x2sq, a_x3sq);
%T_fun = @(x1, x2, x3) a0 + a_x1*x1 + a_x2*x2 + a_x3*x3 + a_x1x2*x1*x2 + a_x1x3*x1*x3 + a_x2x3*x2*x3 + a_x1sq*(x1^2) + a_x2sq*(x2^2) + a_x3sq*(x3^2)
P_rov = sprintf('P = %.4f + %.4fx1 + %.4fx2 + %.4fx3 + %.4fx1x2 + %.4fx1x3 + %.4fx2x3 + %.4fx1^2 + %.4fx2^2 + %.4fx3^2', ...
                    b0, b_x1, b_x2, b_x3, b_x1x2, b_x1x3, b_x2x3, b_x1sq, b_x2sq, b_x3sq);
%P_fun = @(x1, x2, x3) b0 + b_x1*x1 + b_x2*x2 + b_x3*x3 + b_x1x2*x1*x2 + b_x1x3*x1*x3 + b_x2x3*x2*x3 + b_x1sq*(x1^2) + b_x2sq*(x2^2) + b_x3sq*(x3^2)


%% DERIVACE
syms x1 x2 x3 lambda

% Definice funkcí T_fun a P_fun
T_fun = a0 + a_x1*x1 + a_x2*x2 + a_x3*x3 + a_x1x2*x1*x2 + a_x1x3*x1*x3 + a_x2x3*x2*x3 + a_x1sq*(x1^2) + a_x2sq*(x2^2) + a_x3sq*(x3^2);
P_fun = b0 + b_x1*x1 + b_x2*x2 + b_x3*x3 + b_x1x2*x1*x2 + b_x1x3*x1*x3 + b_x2x3*x2*x3 + b_x1sq*(x1^2) + b_x2sq*(x2^2) + b_x3sq*(x3^2);

% Stanovení Lagrangeovy funkce
L=P_fun+lambda*(T_fun-T0);

% Derivace funkcí podle x1, x2, x3
dT_dx1 = diff(T_fun, x1);
dT_dx2 = diff(T_fun, x2);
dT_dx3 = diff(T_fun, x3);

dP_dx1 = diff(P_fun, x1);
dP_dx2 = diff(P_fun, x2);
dP_dx3 = diff(P_fun, x3);

dL_dx1 = diff(L, x1);
dL_dx2 = diff(L, x2);
dL_dx3 = diff(L, x3);
dL_dlambda=diff(L,lambda);


dL2_dx1 = diff(diff(L, 'x1'), 'x1');
dL2_dx2 = diff(diff(L, 'x2'), 'x2');
dL2_dx3 = diff(diff(L, 'x3'), 'x3');
dL2_dlambda=diff(diff(L, 'lambda'), 'lambda');

dL2_dx1dx2 = diff(diff(L, 'x1'), 'x2');
dL2_dx1dx3 = diff(diff(L, 'x1'), 'x3');
dL2_dx2dx3 = diff(diff(L, 'x2'), 'x3');

dL2_dx1dlambda=diff(diff(L, 'x1'), 'lambda');
dL2_dx2dlambda=diff(diff(L, 'x2'), 'lambda');
dL2_dx3dlambda=diff(diff(L, 'x3'), 'lambda');

%% DOSAZENÍ HODNOT
% Počáteční odhady
K_konkretni=-0.7805; 
alfa_konkretni=5.4131;
thetatwist_konkretni=-16.1327;
lambda_konkretni=-4.0756;
X_k=[K_konkretni;alfa_konkretni;thetatwist_konkretni;lambda_konkretni]*0.999;

for i=1:500
%% První derivace
dL_dx1_hodnota=eval(subs(dL_dx1, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL_dx2_hodnota=eval(subs(dL_dx2, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL_dx3_hodnota=eval(subs(dL_dx3, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL_dlambda_hodnota=eval(subs(dL_dlambda, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));


%% Druhé derivace
dL2_dx1_hodnota=eval(subs(dL_dx1, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dx2_hodnota=eval(subs(dL_dx2, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dx3_hodnota=eval(subs(dL_dx3, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dlambda_hodnota=eval(subs(dL2_dlambda, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));

dL2_dx1dx2_hodnota=eval(subs(dL2_dx1dx2, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dx1dx3_hodnota=eval(subs(dL2_dx1dx3, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dx2dx3_hodnota=eval(subs(dL2_dx2dx3, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));

dL2_dx1dlambda_hodnota=eval(subs(dL2_dx1dlambda, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dx2dlambda_hodnota=eval(subs(dL2_dx2dlambda, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));
dL2_dx3dlambda_hodnota=eval(subs(dL2_dx3dlambda, {'x1', 'x2', 'x3','lambda'}, {K_konkretni, alfa_konkretni, thetatwist_konkretni, lambda_konkretni}));

%% Hessián a gradient Lagrange
H=[dL2_dx1_hodnota, dL2_dx1dx2_hodnota, dL2_dx1dx3_hodnota, dL2_dx1dlambda_hodnota;
    dL2_dx1dx2_hodnota, dL2_dx2_hodnota, dL2_dx2dx3_hodnota, dL2_dx2dlambda_hodnota;
    dL2_dx1dx3_hodnota, dL2_dx2dx3_hodnota, dL2_dx3_hodnota, dL2_dx3dlambda_hodnota;
    dL2_dx1dlambda_hodnota, dL2_dx2dlambda_hodnota, dL2_dx3dlambda_hodnota, dL2_dlambda_hodnota];

determinant_H=det(H);

H_inv = inv(H);

vlastni_cisla=eig(H);

Lagrange_g=[dL_dx1_hodnota;dL_dx2_hodnota;dL_dx3_hodnota;dL_dlambda_hodnota];

X_kplus1_hvezd=X_k-H_inv*Lagrange_g;

X_kplus1=X_k+0.1*(X_kplus1_hvezd-X_k);

X_kplus1_vektory(:,i)=X_kplus1;

if X_kplus1(1)>=-0.5 
    X_kplus1(1)=-0.5;
end
if X_kplus1(1)<=-0.92 
    X_kplus1(1)=-0.92;
end

if X_kplus1(2)<=4 
    X_kplus1(2)=4;
end 

if abs(K_konkretni-X_kplus1(1))<=0.0001 && abs(alfa_konkretni-X_kplus1(2))<=0.0001 && abs(thetatwist_konkretni-X_kplus1(3))<=0.0001
    K_konec=X_kplus1(1);
    alfa_konec=X_kplus1(2);
    thetatwist_konec=X_kplus1(3);
    break
else
    K_konkretni = X_kplus1(1);
    alfa_konkretni = X_kplus1(2);
    thetatwist_konkretni = X_kplus1(3);
    lambda_konkretni=X_kplus1(4);
end
iterace=i

if i==500
    disp('dosaženo max. počet iterací')
    break
end

end

% Kontrola regresního modelu
for u=1:length(K)
    T_kontrol(u)=eval(subs(T_fun, {'x1', 'x2', 'x3'}, {K(u), alfa(u), thetatwist(u)}));
    P_kontrol(u)=eval(subs(P_fun, {'x1', 'x2', 'x3'}, {K(u), alfa(u), thetatwist(u)}));
end
fprintf('Tah - kontrolni vypocet\n')
disp(T_kontrol)
rozdil_kontrola_T=T-T_kontrol
fprintf('Vykon - kontrolni vypocet\n')
disp(P_kontrol)
rozdil_kontrola_P=P-P_kontrol

fprintf('Výsledná parametrizace geometrie:\n\n')
fprintf('Zúžení=')
disp(X_kplus1(1))
fprintf('\n')
fprintf('alfa75=')
disp(X_kplus1(2))
fprintf('\n')
fprintf('Zkrouceni=')
disp(X_kplus1(3))