MATLAB: Equality and inequality constraint

Global Optimization ToolboxmogaOptimization Toolbox

Dear All,
Currently, I am doing an integrated biorefinery process of two chemical production.
This is my objective function
p(1) = -(21655.11+61.8*x(1)-906.31*x(2)+9.34*x(2)^2)-(5924.55 -2.16*x(1)…
+ 142.75*x(3)+ 0.42*x(1)*x(3)- 0.0648*x(1)^2-0.911*x(3)^2); %% NPV
p(2) = (278.15+2.30*x(1)-6.38*x(2)+0.017*x(1)*x(2)-0.008*x(1)^2+0.075*x(2)^2)…
-(6.62+0.036*x(1)+0.21*x(3)+0.00016*x(1)*x(3)-0.00012*x(1)^2-0.0012*x(3)^2); %%TDI
p(3) =(9906.56+1660.70*x(1)-345.99*x(2)+0.941*x(1)*x(2)+0.1776*x(1)^2+3.13*x(2)^2)…
+(44362.18+25.12*x(1)^2+15.214466275*x(3)^2); %% GWP
Constraint
function [c,c_eq] = constraints(x)
c = [14.22-0.54*x(2)+0.003*x(1)*x(2)+0.005*x(2)^2-10;…
-40.52-0.022*x(1)+1.03*x(3)+0.0032*x(1)*x(3)-0.0067*x(3)^2-15]; %% This is a constraint of demand
c_eq = [(((14.22-0.54*x(2)+0.003*x(1)*x(2)+0.005*x(2)^2)/x(2))*100)-(((40.52-0.022*x(1)+1.03*x(3)+0.0032*x(1)*x(3)-0.0067*x(3)^2)/x(3))*100)- (0.3*x(1))]; %% Then this is constraint of Total glucose consumption
end
I encounter a difficulty to run my model as it only shows single dot solution. would you like to check whether there is a mistake on my model?

Best Answer

When you have an equality constraint, it is common to be able to get further by solving the equality for one of the variables and substituting that definition for the variable into the other portions of the function.
If you solve ceq for any one of X(1) or X(2) or X(3), you get two solutions -- that is, it is quadratic in each of the variables. We arbitrarily choose to solve ceq for the third variable X(3). We then substitute each of the two solutions in to c, getting one nonlinear inequality constraint for each of the two roots of X(3). We also substitute each of the two roots in to the fitness function, getting one fitness function for each of the two roots of X(3). When then run gamultiobj for each of the two, and then combine the results in a single graph.
Because there are upper and lower bounds on X(3), any pareto front location we find induced by the two situations might or might not meet the bounds criteria when be back-substitute the X(1) and X(2) locations into the derived equations for X(3) and check the result against the bounds. But rather than just discarding the locations that are out of bounds, we can plot them in a different color, to show the pareto front locations that were located and worked, and also the pareto front locations that were located and which did not work. We plot the locations that worked as blue downward triangles for the lower root of X(3) and as blue upward triangles for the upper root of X(3). We plot the locations that fail the X(3) bounds check as red down and upwards triangles for the two roots of X(3).
When you look at the results... you will see only red. This means that the pareto front locations for the first two variables under the equality constraint on X(3), are at places outside the bounds for X(3).
There is no solution for the combined constraints.
You can take this code and modify it to solve for a different variable in hopes that the pareto front gets lucky to find a solution, but you also need to consider the possibility that there are no pareto front locations within the constraints you posted.
I would recommend re-checking the constraint equations.
% clc
% clear  
% close 
%%-----------------------------------------------------------------------%%
% Optimize with gamultiobj
options = optimoptions('gamultiobj','Display','iter',...
          'MaxGeneration',1000,...
          'PopulationSize',50,...
          'CrossoverFraction',0.8,...
          'MigrationFraction',0.01,...
          'PlotFcn',@gaplotpareto);
fitness =@biorefinery;
nvars = 3;
LB = [50 0.67 0.76];
UB = [100 0.77 0.82];
ConsFcn =@constraints; 
[x,fval] = gamultiobj(fitness,nvars,[],[],[],[],LB,UB,ConsFcn,options);
% Plot results
pareto front
figure(1); 
scatter3(fval(:,1),fval(:,2),fval(:,3),'o');
xlabel('NPV ($Million)');
ylabel('FEDI');
zlabel('GWP (kg CO2-eq)');
% view(15,40)
X = sym('X', [1 3]);
[sc, sceq] = ConsFcn(X);
sX3 = solve(sceq, X(3));
scX3a = simplify(subs(sc, X(3), sX3(1)));   %quadratic, select smaller root
scX3b = simplify(subs(sc, X(3), sX3(2)));   %quadratic, select larger root
fscX3a = matlabFunction(scX3a, 'Vars', {X(1:2)});
fscX3b = matlabFunction(scX3b, 'Vars', {X(1:2)});
fitness3a = matlabFunction( subs(fitness(X), X(3), sX3(1)), 'Vars', {X(1:2)});
fitness3b = matlabFunction( subs(fitness(X), X(3), sX3(2)), 'Vars', {X(1:2)});
[x3a, fval3a] = gamultiobj(fitness3a, 2, [], [], [], [], LB(1:2), UB(1:2), @(X) deal(fscX3a(X),[]), options);
[x3b, fval3b] = gamultiobj(fitness3b, 2, [], [], [], [], LB(1:2), UB(1:2), @(X) deal(fscX3b(X),[]), options);
bc3a = double( subs(sX3(1), {X(1), X(2)}, {x3a(:,1), x3a(:,2)}) );
bc3b = double( subs(sX3(2), {X(1), X(2)}, {x3b(:,1), x3b(:,2)}) );
ib3a = bc3a >= LB(3) & bc3a <= UB(3);
ib3b = bc3b >= LB(3) & bc3b <= UB(3);
figure(2)
ps = 20;
ha1 = scatter3(fval3a(ib3a,1),fval3a(ib3a,2),fval3a(ib3a,3), ps, 'b', 'v');
hold on
ha2 = scatter3(fval3a(~ib3a,1),fval3a(~ib3a,2),fval3a(~ib3a,3), ps, 'r', 'v');
hb1 = scatter3(fval3b(ib3b,1),fval3b(ib3b,2),fval3b(ib3b,3), ps, 'b', '^');
hb2 = scatter3(fval3b(~ib3b,1),fval3b(~ib3b,2),fval3b(~ib3b,3), ps, 'r', '^');
hold off
xlabel('NPV ($Million)');
ylabel('FEDI');
zlabel('GWP (kg CO2-eq)');
legend([ha1, ha2, hb1, hb2], {'lower root in bounds', 'lower root out of bounds', 'upper root in bounds', 'upper root out of bounds'})
Related Question