Skip to content
Snippets Groups Projects
Select Git revision
  • 3be59c79733282dc0fe3c29374c24e861aeb8be1
  • master default protected
2 results

Skript.m

  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    Skript.m 4.57 KiB
    %% Matlab file for the Paper 2024:
    % - reproduce the simulations
    % - reproduce optimisation results
    % - reproduce all figures.png (converted to .eps using LibreOfficeDraw)
    % - test average computation time
    
    %% Numerical Simulations
    % figures for SWRO, PRO and the ERD
    % original .m file from: C:>User>bav1839>Documents>MATLAB>Counter-current-solver
    close all;
    [Output1, Output2, time]=fun_1(1,0,'sol',1e6,1e-6),
    %figure(1); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_6.png');
    %figure(2); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_7.png');
    %figure(3); set(gcf,'color','w'); f = gcf; exportgraphics(f,'Figure_8_ERD.png');
    %
    
    %% Average simulation time:
    % manually disable figures first
    T=zeros(1,100);
    for i=1:length(T)
    [Sim_3b_output1, Sim_3b_output2, time]=fun_1(1,.32,'sol',1e6,1e-6);
    T(i)=time;
    end
    mean(T)
    %5.9319
    
    %% Optimisation-1(blue)
    c1=datetime("now");
    x0=[69.9 51.5];
    FWmin=linspace(0.3770, 0.0015, 40);
    for i=2:39
    clc,
    i, 
    close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
    foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
    option_mesh = 1e4; option_BVP = 1e-4; option_data = 1;
    [x_opt, f_opt]=fmincon(@(x)fun_1(x,1,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',1,option_mesh,option_BVP, FWmin(i)),foptions);
    end
    c2=datetime("now"); %von 2 bis 39
    save DATA_T1.mat c1 c2 % 00:02:30
    
    %% Optimisation-2 (red)
    close all;clear all;clc
    c3=datetime("now");
    x0=[69.996767760342830 49.753887785404920];
    FWmin=linspace(0.3862, 0.0066, 40);
    for i=2:39
    clc,
    i,
    close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
    foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-16, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
    option_mesh = 1e4; option_BVP = 1e-4; option_data = 2;
    [x_opt, f_opt]=fmincon(@(x)fun_1(x,2,'SEC',option_mesh,option_BVP),x0,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW',2,option_mesh,option_BVP, FWmin(i)),foptions);
    end
    c4=datetime("now"); %von 2 bis 39
    save DATA_T2.mat c3 c4
    (c4-c3)/38, %00:01:54
    
    %% hybrid paretosearch - old
    close all; clear all; clc
    rng default % For reproducibility
    A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5];
    % initial points:
    load('DATA_Approx.mat');
    X0 = X_approx'; %size(X0)= (93,5)
    %remove infeasible poitns
    z=length(X_approx(1,:));
    for i=1:z
        k=0;
        for j=1:4
        if X_approx(j,z+1-i) > ub(j); k=1;end
        if X_approx(j,z+1-i) < lb(j); k=1;end
        end
        if k==1;X0(z+1-i,:)=[]; end
    end
    %After removing points:
    %size(X0)= (61,5)
    %
    options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000);
    % x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
    [X_paretosearch_3, Y_paretosearch_3] = paretosearch(@(x)fun_1(x,3,'Pareto',1e4,1e-4),5,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); 
    
    %%
    load DATA_Pareto.mat
    scatter(SEC_Pareto_1,FW_Pareto_1);hold on
    scatter(SEC_Pareto_2,FW_Pareto_2)
    
    %% red - Paretosearch
    %
    close all; clear all; clc
    rng default % For reproducibility
    %
    A= [-1 1]; b= 0; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
    %A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5];
    %
    % initial points:
    load DATA_Pareto.mat
    X0 = [X_Pareto_2';X_Pareto_2';X_Pareto_2';X_Pareto_2';X_Pareto_2'];
    %
    options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000);
    c5=datetime("now");
    % x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
    [X, Y] = paretosearch(@(x)fun_1(x,2,'Pareto',1e4,1e-4),2,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); 
    c6=datetime("now");
    scatter(-Y(:,1),-Y(:,2));hold on
    scatter(SEC_Pareto_2,FW_Pareto_2)
    save DATA_T3 c5 c6
    %c6-c5 = 00:56:57
    
    %% blue - Paretosearch
    %
    close all; clear all; clc
    rng default % For reproducibility
    %
    A= [-1 1]; b= 0; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
    %A= []; b= []; Aeq=[]; beq=[]; lb = [30;10;10;1.05;1.5]; ub = [70;20;20;2;1.5];
    %
    % initial points:
    load DATA_Pareto.mat
    X0 = [X_Pareto_1';X_Pareto_1';X_Pareto_1';X_Pareto_1';X_Pareto_1'];
    %
    options = optimoptions('paretosearch','ParetoSetSize',200, 'InitialPoints',X0,'Display','iter', 'MaxFunctionEvaluations',5000);
    c7=datetime("now");
    % x = paretosearch(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
    [X, Y] = paretosearch(@(x)fun_1(x,1,'Pareto',1e4,1e-4),2,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'default'),options); 
    c8=datetime("now");
    scatter(-Y(:,1),-Y(:,2));hold on
    scatter(SEC_Pareto_1,FW_Pareto_1)
    save DATA_T4 c7 c8
    %c8-c7 = 01:05:44