Vectorized the SWRO/PRO-evaluation

Sim_0a_data.m

+function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0a_data(input1, input2)
+%%  Sim_0a_data(input)
+%
+%   Data for Simulation in chapter 3
+%
+%   SWRO-PRO hybrid system I (realistic ERD)
+%
+%   option_data = 0
+%
+%   Input:
+%       input1        -   x
+%       input2        -   x
+
+%% model versions
+version=zeros(1,10);
+% version(1)=0 if co-current, 1 otherwise
+version(2)=1;  % 0 = SWRO beta fixed
+version(3)=0;  % 0 = PRO beta fixed
+version(4)=1;  % 0 = ideal SWRO
+version(5)=1;  % 0 = ideal PRO
+% ERD
+version(6)=3;
+% version(6) = 0 --> only SWRO (no ERD)
+% version(6) = 1 --> only SWRO (with ERD)
+% version(6) = 2 --> only PRO
+% version(6) = 3 --> SWRO-PRO hybrid system (with one ERD)
+% version(6) = 4 --> SWRO-PRO hybrid system (with two ERDs) 
+version(7)=1; % 1 = ICP and ECP for SWRO
+version(8)=1; % 1 = ICP and ECP for PRO
+
+%% Membrane unit properties
+H = 1e-3;           % height of the membrane [m]
+ro_water = 1000;    % mass density of water  [kg/m^3] 
+ro_salt = 2165;     % mass density of salt [kg/m^3] 
+Mw = 18;            % molecular weight of water [kg/kmol]
+Ms = 58.44;         % molecular weight of salt  [kg/kmol]
+Rw = 462;           % gas constant of water  [J/(kg K)] 
+T0 = 297;           % temperature [K] 
+eta = 1.3e-3;       % seawater viscosity [kg/(m s)]
+p_r = 1e5;          % pressure [Pa]=[kg/ms^2]
+rho_r = 1e3;        % density [kg/m^3]
+C_r = 1;            % salt concentration [%]
+
+%% SWRO
+swro_Z=1;               % width of SWRO membrane [m]
+swro_L=4;               % length of SWRO membrane [m]
+if version(6)== 2; swro_L=1; end
+swro_alpha = 5.0815e-9; % SWRO water permeablity co-efficient [s/m]
+swro_R =0.96;           % SWRO salt rejection rate 
+swro_KK = 1e-2;         % SWRO ICP mass transfer coefficient
+swro_KD = 1/swro_KK;    % SWRO ECP draw side mass transfer coefficient
+swro_KF = 1/swro_KK;    % SWRO ECP fresh side mass transfer coefficient
+swro_x_r=swro_L^2;      % x_r=swro_L^2 since x = linspace(0,1,n) (if x = linspace(0,swro_L,n) then x_r=swro_L;)     
+swro_b1 = H/swro_x_r;                           % H/swro_L ratio
+swro_b2 = swro_Z/swro_x_r;                      % Z/swro_L ratio
+J_r = sqrt(H^3/swro_x_r*p_r*rho_r);             % flux [kg/s^2]
+swro_gamma= swro_x_r * p_r *  swro_alpha /J_r;  % SWRO scaling factor - mass balance
+swro_gamma2= J_r^2./(swro_x_r^2 * p_r * rho_r); % SWRO scaling factor - momentum balance
+swro_W_r=J_r*p_r/rho_r;                         % net work [W/m^2]
+sigma=0.999 ;                                   % Rejection coefficient
+swro_beta_fix=4.43e-4/J_r*swro_x_r;             % value for fixed SWRO beta [kg/sm^2]
+ 
+%% PRO
+Z=1;                    % width of the PRO membrane [m]
+L=1.5;                  % length of the PRO membrane [m]      
+alpha = 5.47e-9;        % water permeablity co-efficient [s/m]
+R =0.94;                % salt rejection rate [1]
+KK=7.13e2;              % mass transfer coefficient [sm^2/kg]
+KD=1/KK;                % PRO ECP draw side mass transfer coefficient [sm^2/kg]
+KF=1/KK;                % PRO ECP fresh side mass transfer coefficient [sm^2/kg]
+x_r=L^2;                % x_r=L^2 since x = linspace(0,1,n) (if x = linspace(0,L,n) then x_r=L;)              
+b1 = H/x_r;             % H/L ratio
+b2 = Z/x_r;                           % Z/L ratio
+Q_r = sqrt(H^3/x_r*p_r*rho_r);        % flux [kg/s^2]
+gamma= x_r * p_r * alpha /Q_r ;       % PRO scaling factor - mass balance
+gamma2= Q_r^2./(x_r^2 * p_r * rho_r); % PRO scaling factor - momentum balance
+W_r=Q_r*p_r/rho_r;                    % net work [W/m^2]
+beta_fix =1.71e-4/Q_r*x_r;            % value for fixed PRO beta [kg/sm^2]
+
+%% Sea Water
+cE= 35/983/C_r;                            % salt concentration in seawater
+pE= 1e5/p_r;                               % external pressure
+rho_E=(cE + 1)./(ro_water*cE/ro_salt + 1); % density of incomming seawater
+
+%% SWRO operating conditions
+J_sf_0 = 0;              % salt flux in fresh side at 0 
+J_wf_0 = 0;              % water flux in fresh side at L               
+Pf_L = pE;               % pressure fresh side at L
+
+Pd_0 =  66e5/p_r;        % pressure draw side at 0
+Pd_L =  60e5/p_r;        % pressure draw side at L (not need in the hybrid system)
+
+%% PRO operation conditions 
+Q_sf_0 = 0;              % salt flux in fresh side at 0
+pf_L = pE;               % pressure of fresh side at L
+
+pd_0 = 12e5/p_r;         % pressure draw side at 0
+pd_L = 14.6e5/p_r;       % pressure of fresh side at 0
+pf_0 = 1.3e5/p_r;        % pressure draw side at L 
+
+%% ERD/Turbine/Pump parameters
+T_eff  = .95;               % turbine efficiency
+HP_eff = .9;            	% high pressure pump efficiency
+LP_eff = .95;               % low pressure pump efficiency
+V_m=0.052;                  % Volumetric mixing
+ERD_eff=.96;                % ERD unit pressure efficiency
+ERD_fric=5e-04;             % ERD friction coefficient
+A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
+eta_ERD=0.01;               % leak of high pressure brine
+mix_density=997/rho_r;  	% density of mixture in ERD
+pw=2;                       % water price [$/m^3]
+pe=0.5;                     % electricity price [$/kWh]
+mixer_ERD=1;                % PRO Draw outlet mixer adjustment (only if 2nd ERDs) (mixer_ERD=1 --> all flow to ERD2 no turbine needed)    
+
+%% display figures
+fig=[1,0,0,0,0,0,0,0,1,0,0]; % f(i)=1 --> figure i will be displayed
+%fig=[0,0,0,0,0,0,0,0,0,0,0]; % f(i)=1 --> figure i will be displayed
+
+%% model specific changes:
+% co-current
+if pd_0 > pd_L; version(1)=0; else; version(1)=1; end            
+% no SWRO needed (--> use trivial data)
+if version(6) == 2; version(4)=0; Pd_0=40; Pd_L=35; end
+
+end
\ No newline at end of file

Sim_0_data.m → Sim_0b_data.m

-function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0_data(input1, input2)
-%%  Sim_32_data(input)
+function [H,Z,swro_Z,ro_water,ro_salt,Mw,Ms,Rw,T0,eta,sigma,p_r,rho_r,C_r,swro_L,swro_alpha,swro_R,swro_KK,swro_x_r,swro_b1,swro_b2,J_r,swro_gamma,swro_gamma2,swro_W_r,L,alpha,R,KK,x_r,b1,b2,Q_r,gamma,gamma2,W_r,cE,pE,rho_E,J_sf_0,J_wf_0,Pd_0,Pd_L,Pf_L,Q_sf_0,pd_0,pf_0,pd_L,pf_L,HP_eff,LP_eff,T_eff,V_m,ERD_eff,ERD_fric,A_ERD,eta_ERD,mix_density,pw,pe,swro_beta_fix,beta_fix,mixer_ERD,version,fig,swro_KF,swro_KD,KF,KD]= Sim_0b_data(input1, input2)
+%%  Sim_0b_data(input)
 %
 %   Data for Simulation in chapter 3
 %
 %   SWRO-PRO hybrid system I (realistic ERD)
 %
-%   option_data = 0
+%   option_data = 0.1
 %
 %   Input:
 %       input1        -   x
@@ -104,7 +104,7 @@ HP_eff = .9;            	% high pressure pump efficiency
 LP_eff = .95;               % low pressure pump efficiency
 V_m=0.052;                  % Volumetric mixing
 ERD_eff=.96;                % ERD unit pressure efficiency
-ERD_fric=5e-06;             % ERD friction coefficient
+ERD_fric=5e-04;             % ERD friction coefficient
 A_ERD=H*swro_Z;             % cross sectional area of ERD inflows/outflows
 eta_ERD=0.01;               % leak of high pressure brine
 mix_density=997/rho_r;  	% density of mixture in ERD

Skript.asv

+%% Matlab skript 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
+%
+% Create figures form chapter: Numerical simulations
+% 
+close all;
+[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4),
+%
+%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:
+%
+% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data"
+%
+T=zeros(1,100); 
+for i=1:length(T)
+[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4);
+T(i)=time;
+end
+mean(T) %= 0.9600
+%
+
+%% Optimisation - case2 (with epsilon constraints)
+%
+% Create Paretofront for SWRO+ERD 
+%
+load("DATA_case2.mat")
+% 
+% load("DATA_case2.mat")
+% save DATA_case2.mat X_2
+% 
+% X2  -  set of operating conditions for epsilon constraint method (2x40)
+% Y2  -  system output using X2 (4x40)
+%
+% X2_P  -  set of operating conditions from paretosearch (2x200)
+% Y2_P  -  system output using X2_P (4x200)
+%
+% time2    - time for epsilon constraint method (40points)
+% time2_P  - time for paretosearch solver (200 points)
+%
+
+X2=zeros(2,40);Y2=zeros(4,40);X2_P=zeros(2,200);
+
+
+
+load("DATA_T1.mat")
+c1=datetime("now");
+x0=[69.9 51.5];
+FWmin=linspace(0.3770, 0.0015, 40);
+for i=1:40
+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 - case2 (with paretosearch)
+%
+% Create Paretofront for SWRO+ERD 
+%
+
+
+
+%% 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

Skript.m

@@ -20,21 +20,40 @@ close all;
 
 %% Average simulation time:
 %
-% IMPORTANT: manually comment out line: 116<->117 in "SIM_0_data"
+% calculate average computation time of previous simulation
 %
 T=zeros(1,100); 
 for i=1:length(T)
-[Output1, Output2, time]=fun_1(1,0,'sol',1e4,1e-4);
+[Output1, Output2, time]=fun_1(1,0.1,'sol',1e4,1e-4);
 T(i)=time;
 end
 mean(T) %= 0.9600
 %
 
-%% Optimisation-1(blue)
+%% Optimisation - case2 (with epsilon constraints)
+%
+% Create Paretofront for SWRO+ERD 
+%
+load("DATA_case2.mat")
+% 
+% load("DATA_case2.mat")
+% save DATA_case2.mat X2 Y2 X2_P Y2_P time2 time2_P
+% 
+% X2  -  set of operating conditions for epsilon constraint method (2x40)
+% Y2  -  system output using X2 (4x40)
+%
+% X2_P  -  set of operating conditions from paretosearch (2x200)
+% Y2_P  -  system output using X2_P (4x200)
+%
+% time2    - time for epsilon constraint method (40points)
+% time2_P  - time for paretosearch solver (200 points)
+%
+
+load("DATA_T1.mat")
 c1=datetime("now");
 x0=[69.9 51.5];
 FWmin=linspace(0.3770, 0.0015, 40);
-for i=2:39
+for i=1:40
 clc,
 i, 
 close all; A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
@@ -45,6 +64,13 @@ end
 c2=datetime("now"); %von 2 bis 39
 save DATA_T1.mat c1 c2 % 00:02:30
 
+%% Optimisation - case2 (with paretosearch)
+%
+% Create Paretofront for SWRO+ERD 
+%
+
+
+
 %% Optimisation-2 (red)
 close all;clear all;clc
 c3=datetime("now");