function Aufgabe_2_Optimierung
clc
disp('|')
%if caught error %Stoppt den Run und springt zu der Stelle an welcher der Code
%Probleme verursacht
%% Randbedingungen
disp('Optimierung eines [0/90/+45/-45]s -Laminats hinsichtlich Gewicht/ Schichtdicke')
disp('--------------------------------------------------------------------------')
global k alpha Anz_Lastfaelle n_dach R_pz R_pd R_sz R_sd R_sp Q_12 t_Laminat sigma_12

% Schichtanzahl: [0/90/+45/-45/-45+45/90/0]
k=8;

% Faserrichtung (Referenz zur Laminat-x-Achse)
alpha=[0,90,+45,-45]*pi/180;
for kk=1:k/2
    alpha(k+1-kk)=alpha(kk);
end

% Lastvektor (nur In-plane-Lasten: Kraftflüsse) % Projektnr.:5
Anz_Lastfaelle=3;
n_dach{1}=[860,340,-1045]'; % [nx,ny,nxy]
n_dach{2}=[465,-270,870]';
n_dach{3}=[420,-1005,985]';

% Festigkeitskennwerte
R_pz =1450;        % Zugfestigkeit parallel (in Faserrichtung)
R_pd =-1400;       % Druckfestigkeit parallel (in Faserrichtung)
R_sz =55;          % Zugfestigkeit senkrecht (quer zur Faserrichtung)
R_sd =-1400;       % Druckfestigkeit senkrecht (quer zur Faserrichtung)
R_sp =90;          % Schubfestigkeit

% Materialkonstanten
E_p  =125000;      % E-Modul in Faserrichtung
E_s  =8000;        % E-Modul quer zur Faserrichtung
G_sp =5000;        % Schub-Modul senkrecht/parallel
v_sp =0.3;         % Poisson senkrecht/parallel
v_ps=E_s/E_p*v_sp; % Poisson parallel/senkrecht
% Nachgiebigkeit- bzw. Steifigkeitsmatrix Schicht (2D):
S_12=[1/E_p, -v_ps/E_s, 0;-v_sp/E_p, 1/E_s, 0;0,0,1/G_sp];
Q_12=inv(S_12);

%% Optimierung
% Initialiserungs-Schichtdicken
t_init=[0.1,0.1,0.1];


% lineare Nebenbedingungen
Aeq=[];Beq=[];          % Aeq * t = Beq
A=[];  B=[];            % A   * t ? B
% zulässiges Schichtdicken-Intervall: LB ? t ? UB
LB=[0.01,0.01,0.01];    % untere Grenze: nur positive Werte
UB=[5.00,5.00,5.00];    % obere Grenze
% Funktionsaufruf
[t,t_Laminat]=fmincon(@Zielfunktion,t_init,A,B,Aeq,Beq,LB,UB,@Nebenbedingung);

% Ausgabe
disp('--------------------------------------------------------------------------')
disp('Schichtdicken:')
disp(['    0°-Schicht(t_0) :       ',num2str(t(1)),' mm'])
disp(['   90°-Schicht(t_90):       ',num2str(t(2)),' mm'])
disp([' +-45°-Schicht(t_45):       ',num2str(t(3)),' mm'])
disp(['Laminatdicke (t_Laminat) = 2*(t_0+t_90+2*t_45): ',num2str(t_Laminat),' mm'])

%% Ausgabe Einzelschichtspannungen für optimierte Schichtdicke
disp('--------------------------------------------------------------------------')
disp('Scheibenverzerrungen & Einzelschichtspannungen für optimierte Schichtdicke')
bez = '%2d. &%2d. Schicht (%3d°): % 10.3f   % 10.3f   % 10.3f \n';

for ll=1:Anz_Lastfaelle
    disp('--------------------------------------------------------------------------')
    disp(['Lastfall: ',num2str(ll)])
    disp('Einzelschichtspannungen:')
    disp('  Schicht (Faserwinkel):    Sigma_1      Sigma_2       Tau_21')
    for kk=1:k/2
        fprintf(bez,kk,k+1-kk,alpha(kk)*180/pi,sigma_12{ll}(1,kk),sigma_12{ll}(2,kk),sigma_12{ll}(3,kk));
    end
end

end


function t_Laminat=Zielfunktion(t)
% Schichtdicken des [0/90/+45/-45]s-Laminats; t_+45=t_-45=t_45
t_0 =t(1);
t_90=t(2);
t_45=t(3);
% Gesamtschichtdicke des kompletten Laminats
t_Laminat=2*(t_0+t_90+2*t_45);
end

function [C, Ceq]=Nebenbedingung(t)
global k alpha Anz_Lastfaelle n_dach R_pz R_pd R_sz R_sd R_sp Q_12 t_Laminat sigma_12

%% Schichtdicke
t_0 = t(1);
t_90= t(2);
t_45= t(3);
% Lagenaufbau
t_schicht=[t_0,t_90,t_45,t_45,t_45,t_45,t_90,t_0];

%% Aufbau A-Matrix
% z-Koordianten (Mittelebene: z=0)
z = -t_Laminat/2;
for kk=1:k
    z(kk+1)=z(kk)+t_schicht(kk);
end

% A-Matrix
sum_A=zeros(3,3);
for kk=1:k
    % Transformationen
    Te_xy12{kk}=[(cos(alpha(kk)))^2, (sin(alpha(kk)))^2, 1/2*sin(2*alpha(kk));...
                 (sin(alpha(kk)))^2, (cos(alpha(kk)))^2, -1/2*sin(2*alpha(kk));...
                  -sin(2*alpha(kk)),   sin(2*alpha(kk)), cos(2*alpha(kk))];
    Ts_12xy{kk}=Te_xy12{kk}';
    % Steifigkeitsmatrix Laminat (2D)
    Q_xy{kk}=Ts_12xy{kk}*Q_12*Te_xy12{kk};
    % Elemente der A-Matrix
    for ii=1:3
        for jj=1:3
        A(ii,jj)=sum_A(ii,jj)+Q_xy{kk}(ii,jj)*(z(kk+1)-z(kk));
        sum_A(ii,jj)=A(ii,jj);
        end
    end
end
A_inv=inv(A);

%% Einzelschichtdehnungen und -Spannungen
for ll=1:Anz_Lastfaelle
    % MSV-Verzerrungen epsilon_dach = Scheibenverzerrungen epsilon_0
    epsilon_dach=A_inv*n_dach{ll};
    for kk=1:k/2
        % Einzelschichtdehnungen
        epsilon_12{kk}=Te_xy12{kk}*epsilon_dach;
        % Einzelschichtspannungen
        sigma_12{ll}(:,kk)=Q_12*epsilon_12{kk};
    end
end

%% nichtlineare Ungleichheitsnebenbedingungen
jj=1;   % Iterationsindex
% Schleife über Lastfälle
for ll=1:Anz_Lastfaelle
    % Schleife über eine Laminat-Symmetriehälfte
    for kk=1:k/2                                % k=8: Anzahl Schichten
        C(jj) =   sigma_12{ll}(1,kk)-R_pz;      % 1.NB
        jj=jj+1;
        C(jj) = -(sigma_12{ll}(1,kk)-R_pd);     % 2.NB
        jj=jj+1;
        C(jj) =   sigma_12{ll}(2,kk)-R_sz;      % 3.NB
        jj=jj+1;
        C(jj) = -(sigma_12{ll}(2,kk)-R_sd);     % 4.NB
        jj=jj+1;
        C(jj) = abs(sigma_12{ll}(3,kk))-R_sp;   % 5.NB
        jj=jj+1;
    end
end

%% nichtlineare Gleichheitsnebenbedingungen
Ceq=[];
end