diff --git a/DATA_case1.mat b/DATA_case1.mat
index 004ec10e0d4293447c5ef96ff1b845873931e076..f65feaed2a210fae64321e106a7291e55f5217d1 100644
Binary files a/DATA_case1.mat and b/DATA_case1.mat differ
diff --git a/DATA_case2.mat b/DATA_case2.mat
index e66e8aa3eabf7b54bb1562ff923c4f07abfb300f..e34ea4de08f603b4aa3e5a1111416bbfb74c624e 100644
Binary files a/DATA_case2.mat and b/DATA_case2.mat differ
diff --git a/Skript.m b/Skript.m
index b820493783a50870aabece2c9a31ea3c1a9ba592..15b4176a97a3f608f78c7299bf449939cbc2a7d4 100644
--- a/Skript.m
+++ b/Skript.m
@@ -39,7 +39,7 @@ load("DATA_case1.mat")
% load("DATA_case1.mat")
% save DATA_case1.mat X1 Y1 X1_P Y1_P time1 time1_P
%
-% X1 - set of operating conditions for epsilon constraint method (2x40)
+% X1 - set of operating conditions for epsilon constraint method (2x200)
% Y1 - system output using X2 (5x40)
%
% X1_P - set of operating conditions from paretosearch (2x200)
@@ -52,23 +52,25 @@ startTime=datetime("now");
A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
option_mesh = 1e4; option_BVP = 1e-6; option_data = 1;
foptions = optimoptions('fmincon','Display','off','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
-%foptions = optimoptions('fmincon','Display','iter-detailed','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-3, 'MaxFunEvals',100);
rng default
%
-FWmin=linspace(0.0065, 0.3861, 40);
-X0=[linspace(60,69.99,40);linspace(55,50,40)];
+FWmin=linspace(0.0065, 0.3861, 200);
+X0=[linspace(60,69.99,200);linspace(55,50,200)];
%
-parfor i=1:40
+parfor i=1:200
epsilon=1;
minFW=-FWmin(i);
x0=X0(:,i);
- while epsilon > .9
+ while epsilon > .1
[x_neu, minSEC] = fmincon(@(x)fun_1(x,option_data,'SEC',option_mesh,option_BVP),x0 ,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW' ,option_data,option_mesh,option_BVP,-minFW ),foptions);
[x0, minFW] = fmincon(@(x)fun_1(x,option_data,'FW' ,option_mesh,option_BVP),x_neu,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'SEC',option_data,option_mesh,option_BVP,-minSEC),foptions);
- epsilon = max(norm(x0-x_neu)),
+ y_neu = fun_1(x_neu,option_data,'Pareto',option_mesh,option_BVP),
+ y_0 = fun_1(x0,option_data,'Pareto',option_mesh,option_BVP),
+ epsilon = min(norm(x0-x_neu), norm(y_neu-y_0));
end
- X1(:,i)=x0; %x_neu;
+ X1(:,i)=x0;
Y1(:,i)=fun_1([X1(:,i)],option_data,'sol',option_mesh,option_BVP);
+ i,
end
endTime=datetime("now");
time1 = endTime - startTime;
@@ -79,7 +81,7 @@ scatter(Y1(1,:),Y1(2,:),'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','r'); hold o
grid on; xlim([-5.501 -.5]); ylim([0 0.48]);view(2);
legend('Single SWRO unit','SWRO with ERD','Location', 'Northeast');ylabel('FW [m^3/h]','FontSize',16);xlabel('SEC_{net} [kWh/m^3]','FontSize',16);
-% Optimisation - case2 (epsilon constraint method)
+%% Optimisation - case2 (epsilon constraint method)
%
% Create Paretofront for SWRO+ERD
%
@@ -99,24 +101,27 @@ load("DATA_case2.mat")
%
startTime=datetime("now");
A= [-1 1]; b= -.1; Aeq=[]; beq=[]; lb = [30;30]; ub = [70;70];
-option_mesh = 1e4; option_BVP = 1e-6; option_data = 1;
+option_mesh = 1e4; option_BVP = 1e-6; option_data = 2;
foptions = optimoptions('fmincon','Display','off','Algorithm','interior-point', 'StepTolerance', 1e-12, 'OptimalityTolerance',1e-4, 'MaxFunEvals',100);
rng default
%
-FWmin=linspace(0.0015, 0.3769, 40);
-X0=[linspace(37,70,40);linspace(36,50,40)];
+FWmin=linspace(0.0015, 0.3769, 16);
+X0=[linspace(37,70,16);linspace(36,50,16)];
%
-parfor i=1:40
+parfor i=1:16
epsilon=1;
minFW=-FWmin(i);
x0=X0(:,i);
- while epsilon > .9
+ while epsilon > .1
[x_neu, minSEC] = fmincon(@(x)fun_1(x,option_data,'SEC',option_mesh,option_BVP),x0 ,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'FW' ,option_data,option_mesh,option_BVP,-minFW ),foptions);
[x0, minFW] = fmincon(@(x)fun_1(x,option_data,'FW' ,option_mesh,option_BVP),x_neu,A,b,Aeq,beq,lb,ub,@(x)nonlcon(x,'SEC',option_data,option_mesh,option_BVP,-minSEC),foptions);
- epsilon = max(norm(x0-x_neu)),
+ y_neu = fun_1(x_neu,option_data,'Pareto',option_mesh,option_BVP),
+ y_0 = fun_1(x0,option_data,'Pareto',option_mesh,option_BVP),
+ epsilon = min(norm(x0-x_neu), norm(y_neu-y_0));
end
- X2(:,i)=x0; %x_neu;
- Y2(:,i)=fun_1([X1(:,i)],option_data,'sol',option_mesh,option_BVP);
+ X2(:,i)=x0;
+ Y2(:,i)=fun_1([X2(:,i)],option_data,'sol',option_mesh,option_BVP);
+ i,
end
endTime=datetime("now");
time2 = endTime - startTime;
@@ -126,8 +131,16 @@ scatter(Y2(1,:),Y2(2,:),'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','b'); hold o
grid on; xlim([-5.501 -.5]); ylim([0 0.48]);view(2);
legend('Single SWRO unit','SWRO with ERD','Location', 'Northeast');ylabel('FW [m^3/h]','FontSize',16);xlabel('SEC_{net} [kWh/m^3]','FontSize',16);
+%% test
-
+X0=[linspace(37,70,16);linspace(36,50,16)];
+A=zeros(1,16);B=A;
+parfor i=1:16
+ Y0=fun_1([X0(:,i)],option_data,'sol',option_mesh,option_BVP);
+ A(i)=Y0(1);
+ B(i)=Y0(2);
+end
+scatter(A,B,'MarkerEdgeColor',[0 0 0],'MarkerFaceColor','g'); hold on
diff --git a/fun_1.m b/fun_1.m
index 8fdab10d8129818059b53266e69a9ae28e0e0d7d..7179491c85a441412c1c141effa8bbaef195d21c 100644
--- a/fun_1.m
+++ b/fun_1.m
@@ -157,20 +157,20 @@ SWRO_Recovery =(J_f(end)./J_d(1))*100; % SWRO freshwater recovery [%]
PRO_Recovery = (1- Q_f(end)./Q_f(1))*100; % PRO recovery rate [%]
% warning flags:
if FW < 0
- fprintf(2,' \nERROR: Waring: Freshwater production is negative! \n');
+ %fprintf(2,' \nERROR: Waring: Freshwater production is negative! \n');
end
if SWRO_Recovery > 100
- fprintf(2,' \nERROR: Waring: SWRO recovery is greater then 100 %% \n');
+ %fprintf(2,' \nERROR: Waring: SWRO recovery is greater then 100 %% \n');
end
if SWRO_Recovery < 0
- fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n');
+ %fprintf(2,' \nERROR: Waring: SWRO recovery is negative! \n');
end
if version(6)>1
if PRO_Recovery > 100
- fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
+ %fprintf(2,' \nERROR: Waring: PRO recovery is greater then 100 %% \n');
end
if PRO_Recovery < 0
- fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
+ %fprintf(2,' \nERROR: Waring: PRO recovery is negative! \n');
end
end
%% Output of System
@@ -199,7 +199,7 @@ W_p4 = 1/LP_eff * (p_f(1)-pE)*(Q_f(1)*Z)./local_ro_f(1); W_p4 = W_p4*W_r;
if pERD < 0; pERD=pMax; end
W_p3 = 1/HP_eff * (P_d(1)-pERD)*(J_ERD*swro_Z)/rho_ERD; W_p3=W_p3*swro_W_r; W_p4=0; W_t=0;
if W_p3 < 0
- fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n');
+ %fprintf(2,' \nERROR: Waring: Pump 3 generates generates power! \n');
sol.stats.maxerr =1000;
end
end