% M-File zum Starten der Modellsimulation 

clear all
close all
clc

% Unterordner einbinden:
path(genpath(pwd),path);
path(genpath('..\..\easyDoE'),path);
path(genpath('..\..\SubFCT'),path);
path(genpath('..\..\Funktionen'),path);
path(genpath('..\Simulationsdaten'),path);

% Laden der MAT-Files für die Streckenmodelle und Ableitungsmodelle:
load('modLOLIMOT_3In_CA50');
load('modLOLIMOT_3In_IMEP');
load('modLOLIMOT_3In_DPMAX');
load('modLOLIMOT_3In_HC');
load('modLOLIMOT_3In_NOX');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Manuelle Einstellungen
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Simulationsschrittweite
stepSize = 0.06;  % Bei 2000 Umdrehungen in der Minute dauert ein Arbeitsspiel (2 Umdrehungen) 0.06sec an!

% Dauer für einen Sprung, etc.! (Innerhalb der Blöcke: CA50_Vari,
% IMEP_Vari und NOX_Vari zu finden)
Zeiteinheit = 150; % Je kleiner die Zeiteinheit gewählt wird, desto weniger Zeit steht für die Sprünge, Rampen, Sinus-Verläufe bereit!
                   % Mit der Zeiteinheit von 150 stehen stets 9sec. je Sprung, Rampe, Sinus bereit, bevor es zur nächsten Veränderung kommt!

% Simulationsdauer vordefinieren
simTime     = 50; % simTime berechnet sich wie folgt: 95*Zeiteinheit*stepSize+T_Sprung!
                   % Zudem sollte simTime durch stepSize teilbar sein!

% PT1-Zeitkonstante 
T_PT1    = 2; % für AGRstat (Stellgröße) --> AGRdyn (Zustandsgröße)

% Anzahl Regelgrößen:
Regel = 3;

% Anzahl Stellgrößen:
Stell = 3;

% Anzahl der über der Zeit betrachteten Größen (hier: Regel-,
% Soll-Stell-, Ist-Stellgrößen und Sollwerte):
ywuU   = 4;

% Stellgrößen entweder für offenen Regelkreis (ohne MP-Regler) oder für Betriebspunkt vor dem Sprung:
SOIstart = 20; % °KWvZOT
FMIstart = 15; % mm^3/Hub
AGRstart = 30; % % 

% Betriebspunkt nach dem Sprung (Wird hier nicht gebraucht):
CA50 = 6;
IMEP = 4;
NOx  = 100;

% "Sollwertsprung":
T_Sprung = 12; % T_Sprung minimal gleich der Einschwingzeit

% Constraints/Beschränkungen:
% Stellgrößen/Zustandsgrößen:
SOImin   = 10; % °KWvZOT
SOImax   = 50;
FMImin   = 10; % mm^3/Hub
FMImax   = 20;
AGRmin   = 10; % %
AGRmax   = 50;

%Regelgrößen: % Regelgrößen-Beschränkungen bislang inaktiv!
CA50min  = -10;
CA50max  = 10;
IMEPmin  = 1;
IMEPmax  = 6;
NOxmin   = 0;
NOxmax   = 200;

% Stellhorizonte (Hp, Hw, Hu) müssen innerhalb des MP-Reglers definiert werden.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Auswahl Optimierungsdatei:

Select_Opt = [23];

% h = 1;
for h = 1:length(Select_Opt)
    
    switch Select_Opt(h)
        case 1 
              T_exp = 0; % Selbst wenn keine Referenztrajektorien verwendet werden, muss der Parameter übergeben werden! Deswegen trotzdem DEFINIEREN!
              Texp_CA50 = 0;
              Texp_IMEP = 0;
              Texp_NOX  = 0;
              SimModell = ('');  
%             load 
%             RGkonstW1 = RGkonstW(1,8:13); % Wichtungen sind konstant für alle 350 Betriebspunkte.
%             load 
%             RGkonstW2 = RGkonstW(1,8:13);
%             load 
%             RGkonstW3 = RGkonstW(1,8:13);
        case 2
              T_exp = 0; 
              Texp_CA50 = 0;
              Texp_IMEP = 0;
              Texp_NOX  = 0;
              SimModell = ('');
              load 2_RG_OptKonst_RQ_Case2_QP_Hp60_Hu3_Fall1_NOXmax200
              RGkonstW1 = RGkonstW(1,8:13);
              load 2_RG_OptKonst_RQ_Case2_QP_Hp60_Hu3_Fall2_NOXmax200
              RGkonstW2 = RGkonstW(1,8:13);
              load 2_RG_OptKonst_RQ_Case2_QP_Hp60_Hu3_Fall3_NOXmax200
              RGkonstW3 = RGkonstW(1,8:13);
%             load 2_RG_OptKonst_RQ_Case2_QP_Hp60_Hu3_Fall3_NOXmax450
%             RGkonstW3 = RGkonstW(1,8:13);  
        case 3
            T_exp = 0; 
            Texp_CA50 = 0;
            Texp_IMEP = 0;
            Texp_NOX  = 0;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
            load 3_RG_OptKonst_RQ_Case1_QP_Hp60_Hu3_Fall3_NOXmax200
            RGkonstW3 = RGkonstW(1,8:13);
        case 4
            T_exp = 0; 
            Texp_CA50 = 0;
            Texp_IMEP = 0;
            Texp_NOX  = 0;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
            load 4_RG_OptKonst_RQ_Case3_QP_Ref0_5_Hp60_Hu3_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
        case 5
            T_exp = 0.5;
            Texp_CA50 = 0.5;
            Texp_IMEP = 0.5;
            Texp_NOX  = 0.5;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
            load 5_RG_OptKonst_RQ_Case3_QP_Ref0_2_Hp60_Hu3_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
        case 6
            T_exp = 0.2;
            Texp_CA50 = 0.2;
            Texp_IMEP = 0.2;
            Texp_NOX  = 0.2;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
            load 6_RG_OptKonst_RQ_Case3_WIMEP3_QP_Ref0_2_Hp60_Hu3_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
        case 7
            T_exp = 0;
            Texp_CA50 = 0;
            Texp_IMEP = 0;
            Texp_NOX  = 0;
            SimModell = ('');
            load 7_RG_OptKonst_RQ_Case2_QP_NormMPC_Hp60_Hu3_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load 7_RG_OptKonst_RQ_Case2_QP_NormMPC_Hp60_Hu3_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load 7_RG_OptKonst_RQ_Case2_QP_NormMPC_Hp60_Hu3_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
        case 8
            T_exp = 0.5;
            Texp_CA50 = 0.5;
            Texp_IMEP = 0.5;
            Texp_NOX  = 0.5;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 9
            T_exp = 1;
            Texp_CA50 = 1;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 10
            T_exp = 0.5;
            Texp_CA50 = 0.5;
            Texp_IMEP = 0.5;
            Texp_NOX  = 0.5;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 11
            T_exp = 1;
            Texp_CA50 = 1;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 12
            T_exp = 0;
            Texp_CA50 = 0;
            Texp_IMEP = 0;
            Texp_NOX  = 0;
            SimModell = ('');
%           load
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 13
            T_exp = 0.2;
            Texp_CA50 = 0.2;
            Texp_IMEP = 0.2;
            Texp_NOX  = 0.2;
            SimModell = ('');
            load 13_RG_OptKonst_RQ_Case2_QP_NormREFmpc0_2_Hp90_Hu3_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load 13_RG_OptKonst_RQ_Case2_QP_NormREFmpc0_2_Hp90_Hu3_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load 13_RG_OptKonst_RQ_Case2_QP_NormREFmpc0_2_Hp90_Hu3_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
        case 14
            T_exp = 0.5;
            Texp_CA50 = 0.5;
            Texp_IMEP = 0.5;
            Texp_NOX  = 0.5;
            SimModell = ('');    
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 15
            T_exp = 1;
            Texp_CA50 = 1;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('');
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
       case 16
            T_exp = 0.2;
            Texp_CA50 = 0.2;
            Texp_IMEP = 0.2;
            Texp_NOX  = 1;
            SimModell = ('');
            load 16_RG_OptKonst_RQ_Case3_QP_NormREFmpcGuete0_5_1_0_05_Hp20_Hu10_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load 16_RG_OptKonst_RQ_Case3_QP_NormREFmpcGuete0_5_1_0_05_Hp20_Hu10_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load 16_RG_OptKonst_RQ_Case3_QP_NormREFmpcGuete0_5_1_0_05_Hp20_Hu10_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
       case 17
            T_exp = 0.2;
            Texp_CA50 = 0.2;
            Texp_IMEP = 0.2;
            Texp_NOX  = 1;
            SimModell = ('');   
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
            
% Messungen nach erfolgreicher BUG-Suche:

       case 19
            T_exp = 1;
            Texp_CA50 = 1;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('OptFile_19_3x3_n3_LIN_Nox_QP_Norm_Ref_SprngRmpSns');
            load OptFile_19_RG_OptKonst_RQ_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load OptFile_19_RG_OptKonst_RQ_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load OptFile_19_RG_OptKonst_RQ_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
       case 20
            T_exp = 1;
            Texp_CA50 = 1;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('OptFile_20_23_3x3_n3_LIN_Nox_QP_Norm_Ref_SprngRmpSns'); 
            load OptFile_20_RG_OptKonst_RQ_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load OptFile_20_RG_OptKonst_RQ_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load OptFile_20_RG_OptKonst_RQ_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
       case 21
            T_exp = 1;
            Texp_CA50 = 0.5;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('OptFile_20_23_3x3_n3_LIN_Nox_QP_Norm_Ref_SprngRmpSns');  
            load OptFile_21_RG_OptKonst_RQ_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load OptFile_21_RG_OptKonst_RQ_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load OptFile_21_RG_OptKonst_RQ_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
        case 22
            T_exp = 0.5;
            Texp_CA50 = 0.5;
            Texp_IMEP = 1;
            Texp_NOX  = 1;
            SimModell = ('OptFile_20_23_3x3_n3_LIN_Nox_QP_Norm_Ref_SprngRmpSns');
            load OptFile_22_RG_OptKonst_RQ_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load OptFile_22_RG_OptKonst_RQ_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load OptFile_22_RG_OptKonst_RQ_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
       case 23
            T_exp = 0.05;
            Texp_CA50 = 0.05;
            Texp_IMEP = 0.05;
            Texp_NOX  = 1;
            SimModell = ('OptFile_20_23_3x3_n3_LIN_Nox_QP_Norm_Ref_SprngRmpSns'); 
            load OptFile_23_RG_OptKonst_RQ_Fall1
            RGkonstW1 = RGkonstW(1,8:13);
            load OptFile_23_RG_OptKonst_RQ_Fall2
            RGkonstW2 = RGkonstW(1,8:13);
            load OptFile_23_RG_OptKonst_RQ_Fall3
            RGkonstW3 = RGkonstW(1,8:13);
       case 24
            T_exp = 0.05;
            Texp_CA50 = 0.05;
            Texp_IMEP = 0.05;
            Texp_NOX  = 1;
            SimModell = ('');  
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);
        case 25
            T_exp = 0.05;
            Texp_CA50 = 0.05;
            Texp_IMEP = 0.05;
            Texp_NOX  = 1;
            SimModell = ('');      
%           load 
%           RGkonstW1 = RGkonstW(1,8:13);
%           load 
%           RGkonstW2 = RGkonstW(1,8:13);
%           load 
%           RGkonstW3 = RGkonstW(1,8:13);

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% Auswahl konstante Wichtungen:
    Select_Typ_KonstW = [1];
    % RGkonstW1 = [];
    RGkonstW2 = [];
    RGkonstW3 = [];
        % Alle 3 Fälle ausgesucht 'Select_Typ_KonstW = [1 2 3];' --> RGkonstW1=[], RGkonstW2=[] & RGkonstW3=[] inaktiv!
        % Nur Fall 1 & 3 ausgesucht 'Select_Typ_KonstW = [1 3];' --> Nur RGkonstW1=[] & RGkonstW3=[] inaktiv ...
        % ... dadurch hat RGkonstW_ins nur 2 Zeilen und die m-Schleife 2
        % Durchgänge und mit Select_Typ_KonstW(m) in eval (siehe Speicherung)
        % die richtige Verknüpfung!
    RGkonstW_ins = [RGkonstW1;RGkonstW2;RGkonstW3];     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for m = 1:length(Select_Typ_KonstW)    

    R1 = RGkonstW_ins(m,1);
    R2 = RGkonstW_ins(m,2);
    R3 = RGkonstW_ins(m,3);
    Q1 = RGkonstW_ins(m,4);
    Q2 = RGkonstW_ins(m,5);
    Q3 = RGkonstW_ins(m,6);
    
AbtastI = (simTime/stepSize)+1;
Sum_Eintraege = AbtastI*Regel; % Regel-, Soll-Stell-, Ist-Stellgrößen und Sollwerte werden während 
                                   % der Simulation alternierend untereinander abgespeichert!
Mat_Ges_SRS = zeros(9+Sum_Eintraege,ywuU); % Die 9 ist Platzhalter für den Inhalt des Vektors 'Einstellungen'!
length_Mat_Ges_SRS = length(Mat_Ges_SRS);

Einstellungen = [0;0;0;R1;R2;R3;Q1;Q2;Q3];
length_Ein = length(Einstellungen);

simOut=sim(SimModell); 
 
Regelgroessen=RegelGr.signals.values;
Sollwerte=SollWerte.signals.values;
SOLLStellgroessen=SollStellGr.signals.values;
ISTStellgroessen=IstStellGr.signals.values;

for r=1:length(Regelgroessen)

    Regel_altern(((r-1)*Regel+1):((r-1)*Regel+Regel),1) = Regelgroessen(1:Regel,1,r);
    SOLLStell_altern(((r-1)*Stell+1):((r-1)*Stell+Stell),1) = SOLLStellgroessen(1:Stell,1,r);
    Soll_altern(((r-1)*Regel+1):((r-1)*Regel+Regel),1)  = Sollwerte(r,:);
    ISTStell_altern(((r-1)*Regel+1):((r-1)*Regel+Regel),1)  = ISTStellgroessen(r,:);
    
end

Mat_Ges_SRS(1:length_Ein,1) = Einstellungen; 
Mat_Ges_SRS(length_Ein+1:length_Mat_Ges_SRS,1) = Regel_altern;
Mat_Ges_SRS(length_Ein+1:length_Mat_Ges_SRS,2) = Soll_altern;
Mat_Ges_SRS(length_Ein+1:length_Mat_Ges_SRS,3) = SOLLStell_altern;
Mat_Ges_SRS(length_Ein+1:length_Mat_Ges_SRS,4) = ISTStell_altern;

eval(['save OptFile' num2str(Select_Opt(h)) '_SprngRmpSns_Fall' num2str(Select_Typ_KonstW(m)) '_x' num2str(Zeiteinheit) '.mat Mat_Ges_SRS']) 

end

end