Fmincon finds different solutions for optimization problem in dependance of initial values

Hello,

I have got a non-linear optimization problem and want to use fmincon to solve it. Even if I guess "good" inital values I do not get the right solution. I just find local solutions depending on the inital values. The solution I want to find is the following:

m_CO: 23.0468 kg

m_CO2: 95.7836 kg

m_H2: 0.3538 kg

m_H2O: 14.8158 kg

calcualted minimum for G: -1.242 GJ

Could anyone help me to find the error in that code?

% considering chemical reactions
% system of reaction basis 1 solved with fmincon for non-linear function
% with eqality and inequality constraint
% components: H2O, CO2, CO and H2, p=1bar and T=1300 K
clear all
clc

% reaction basis = number of species - number of elements
R = 8.3145; % in J/mol/K
T = 1300;   % in K
p = 1e5;    % in Pa
p0 = 1e5;   % in Pa

% at the moment works only for T=1300 K
[CO2, H2O, CO, H2] = read_component_data(T);

% molar mass of "elements" in kg/kmol
C.M = 12;
O.M = 16;
O2.M = 32;

CO2.m = CO2.n0 * CO2.M; 
CO.m = CO.n0 * CO.M;
H2O.m = H2O.n0 * H2O.M; 
H2.m = H2.n0 * H2.M; 

% initial values in kg (for m0) and kg/kmol (for M0_tot)
CO2.m0 = 90;
CO.m0 = 23;
H2.m0 = 35;
H2O.m0 = 15;
M0_tot = 30;

% overall mass
m = CO2.n0 * CO2.M + H2.n0 * H2.M;

% array m_ of initial values
%m_ = [CO2.m0, CO.m0, H2O.m0, H2.m0, M0_tot];
m_ = [96 23 15 0.3 30 ];
[c, ceq] = constraints(m_, C, O, O2, CO, CO2, H2O, H2, m);

% objective function G_min in function
G_min = @(m_)objective(m_, CO, CO2, H2O, H2, m, T, p, p0, R);

A = []; b = []; Aeq = []; beq = []; lb = []; ub = [];

% solver fmincon for non-linear function with eqality and inequality
% constraints 
nonlincon = @(m_) constraints(m_, C, O, O2, CO, CO2, H2O, H2, m);
[x, fval, ef, output, lambda] = fmincon(G_min, m_, A, b, Aeq, beq, lb, ub, nonlincon);

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

fval

fval = -1.2420e+09

m_CO2_sol = x(1);
m_CO_sol = x(2);
m_H2O_sol = x(3);
m_H2_sol = x(4);
M_tot_sol = x(5);

disp(['solver message: ', num2str(output.message)])

solver message: Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>

Optimization completed: The relative first-order optimality measure, 7.602634e-07,
is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint
violation, 1.386979e-11, is less than options.ConstraintTolerance = 1.000000e-06.

disp(['No of iterations: ', num2str(output.iterations)])

No of iterations: 6

disp(['Solution for m(CO2): ', num2str(m_CO2_sol), ' kg, solution for m(CO): ', num2str(m_CO_sol), ' kg, solution for m(H2O): ', num2str(m_H2O_sol), ' kg, solution for m_H2: ', num2str(m_H2_sol), ' kg'])

Solution for m(CO2): 95.7837 kg, solution for m(CO): 23.0468 kg, solution for m(H2O): 14.8158 kg, solution for m_H2: 0.3538 kg

function G_min = objective(m_, CO, CO2, H2O, H2, m, T, p, p0, R)
% unknown masses m_CO2, m_CO, m_H2O, m_H2 and molar mass M_tot of mixture
% inital values in array m_ in kg and kg/kmol

% initial values
m_CO2 = m_(1);
m_CO = m_(2); 
m_H2O = m_(3);
m_H2 = m_(4);
M_tot = m_(5);

% minimizing Gibbs free enthalpy
G(1) = m_CO2 * (CO2.Delta_f_G0/CO2.M + R * T/CO2.M * log((m_CO2 * M_tot/(m * CO2.M)) * p/p0))*10^3;
G(4) = m_H2 * (H2.Delta_f_G0/H2.M + R * T/H2.M * log((m_H2 * M_tot/(m * H2.M)) * p/p0))*10^3; 
G(2) = m_CO * (CO.Delta_f_G0/CO.M + R * T/CO.M * log((m_CO * M_tot/(m * CO.M)) * p/p0))*10^3; 
G(3) = m_H2O * (H2O.Delta_f_G0/H2O.M + R * T/H2O.M * log((m_H2O * M_tot/(m * H2O.M)) * p/p0))*10^3; 

% objective function Gibbs enthalpy G_min of the system
G_min = G(1) + G(2) + G(3) + G(4);
end

function [c, ceq] = constraints(m_, C, O, O2, CO, CO2, H2O, H2, m)
% inequality constraints c and equality constraints ceq 

% initial values
m_CO2 = m_(1);
m_CO = m_(2); 
m_H2O = m_(3);
m_H2 = m_(4);
M_tot = m_(5);

% equality constraints: C balance, O balance, H balance, molar mass of
% mixture
ceq(1) = C.M/CO2.M * m_CO2 + C.M/CO.M * m_CO - C.M/CO2.M * CO2.n0 * CO2.M;
ceq(2) = O2.M/CO2.M * m_CO2 + O.M/CO.M * m_CO + O.M/H2O.M * m_H2O - O2.M/CO2.M * CO2.n0 * CO2.M;
ceq(3) = m_H2 + H2.M/H2O.M * m_H2O - H2.n0 * H2.M;
ceq(4) = 1/M_tot - m_CO2/(m * CO2.M) - m_H2/(m * H2.M) - m_CO/(m * CO.M) - m_H2O/(m * H2O.M);

% inequality constraints: non-negative mass and molar mass
c(1) = 0 - m_CO2;
c(2) = 0 - m_CO;
c(3) = 0 - m_H2;
c(4) = 0 - m_H2;
c(5) = 0 - M_tot;
end

function [CO2, H2O, CO, H2] = read_component_data(T)
% read component data from excel sheet in same folder to handle Delta_f_G0
% at different temperatures
% n0 = 2, in kmol

    table = readtable('comp_data.xlsx');
    % 1300 K
    if T == 1300
        n = 4;
    end
    
    CO2.n0 = table2array(table(2,"CO2"));
    H2.n0 = table2array(table(2,"H2"));
    CO.n0 = table2array(table(2,"CO"));
    H2O.n0 = table2array(table(2,"H2O"));
    
    % molar mass of species in table row 1 in kg/kmol
    CO2.M = table2array(table(1,"CO2"));
    % CH4.M = 16;
    H2O.M = table2array(table(1,"H2O"));
    CO.M = table2array(table(1,"CO"));
    H2.M = table2array(table(1,"H2"));
    
    % TODO: element structures need to be returned
    % molar mass of elements in kg/kmol
    C.M = 12;
    O.M = 16;
    O2.M = 32;
    
    H2O.Delta_f_G0 = table2array(table(n,"H2O"));
    CO.Delta_f_G0 = table2array(table(n,"CO"));
    CO2.Delta_f_G0 = table2array(table(n,"CO2"));
    H2.Delta_f_G0 = table2array(table(n,"H2"));
end

Answer 1

[Torsten]

It's not unusual that the solution for optimization problems depends on the initial guesses for the optimization variables.

If given initial values near the "true" solution, your code is able to reproduce this solution (see above).

You could try "MultiStart" if you can get this solution also if the initial values are further away.

Answer 2

[Matt J]

You might find it easier to reach a good solution if you simplified the implementation of your constraints. In particular, your equality constraints are nonlinear as far as I can tell only because of the term 1/M_tot. If you were to just redefine the fifth unknown as m_(5)=1/M_tot, then these would become linear equality constraints, which you would implement using the Aeq,beq matrix inputs.

% equality constraints: C balance, O balance, H balance, molar mass of
% mixture
ceq(1) = C.M/CO2.M * m_CO2 + C.M/CO.M * m_CO - C.M/CO2.M * CO2.n0 * CO2.M;
ceq(2) = O2.M/CO2.M * m_CO2 + O.M/CO.M * m_CO + O.M/H2O.M * m_H2O - O2.M/CO2.M * CO2.n0 * CO2.M;
ceq(3) = m_H2 + H2.M/H2O.M * m_H2O - H2.n0 * H2.M;
ceq(4) = 1/M_tot - m_CO2/(m * CO2.M) - m_H2/(m * H2.M) - m_CO/(m * CO.M) - m_H2O/(m * H2O.M);

Similarly, positivity bounds like the ones below should not be implemented as nonlinear inequalities. You should just set lb=zeros(5,1).

% inequality constraints: non-negative mass and molar mass
c(1) = 0 - m_CO2;
c(2) = 0 - m_CO;
c(3) = 0 - m_H2;
c(4) = 0 - m_H2;
c(5) = 0 - M_tot;

[Caroline Otto]

I tried it with linear equality constraints and lower bounds, but I don not get the right result neither. Here is my code for the new main:

clear
clc
format long
R = 8.3145; % gas constant in J/mol/K
T = 1300;   % temperature in K
p = 1e5;    % system pressure in Pa
p0 = 1e5;   % reference pressure in Pa
% at the moment T=1300 K, TODO: varying temperature, to be continued in Excel sheet
[CO2, CO, H2O, H2] = read_component_data(T);
% molar mass of "elements" (for constraints) in kg/kmol
C.M = 12;
O.M = 16;
O2.M = 32;
% mass of species in kg calculated with beginning amount of substance n0
% in kmol
CO2.m = CO2.n0 * CO2.M; 
CO.m = CO.n0 * CO.M;
H2O.m = H2O.n0 * H2O.M; 
H2.m = H2.n0 * H2.M; 
% initial values in kg (for m_init) and kg/kmol (for M_tot_init) for optimization
% problem
CO2.m_init = 96;
CO.m_init = 23;
H2O.m_init = 35;
H2.m_init = 15;
M_tot_init = 30;
m_init = CO2.m_init + CO.m_init + H2O.m_init + H2.m_init;
% overall mass in kg
m = CO2.m + CO.m + H2O.m + H2.m;
% array m_ of initial values for unknowns
m_ = [CO2.m_init, CO.m_init, H2O.m_init, H2.m_init, 1/M_tot_init];
% inequality constraints c and equality constraints ceq
% [c, ceq] = constraints(m_, CO2, CO, H2O, H2, O, O2, C);
% objective function G_min
G = objective(m_, CO2, CO, H2O, H2, m, T, p, p0, R);
G_min = @(m_)G;
% Aeq = [0.63636 1 0 0 0;
%     0.4091 0 1 0 0;
%     -0.0454 0 0 1 0;
%     -1/(m * CO2.M) -1/(m * CO.M) -1/(m * H2O.M) -1/(m * H2.M) 1;
%     1 1 1 1 0];
% 
% beq = [84; 54; -4; 0; m];
Aeq = [12/44 12/28 0 0 0;
    32/44 16/28 16/18 0 0;
    0 0 2/18 1 0;
    -1/(m * CO2.M) -1/(m * CO.M) -1/(m * H2O.M) -1/(m * H2.M) 1;
    1 1 1 1 0];
beq = [12/44*CO2.m; 32/44*CO.m; H2.m; 0; m];
A = []; 
b = []; 
% non-negtive mass and molar mass
lb = zeros(5,1); ub = [];
% solver fmincon for non-linear function with eqality and inequality
% constraints 
% nonlincon = @(m_) constraints(m_, CO2, CO, H2O, H2, O, O2, CO2);
options = optimoptions("fmincon");
options.Display = "iter";
[x, fval, ef, output, lambda] = fmincon(G_min, m_, A, b, Aeq, beq, lb, ub, [], options);
CO2.m_sol = x(1);
CO.m_sol = x(2);
H2O.m_sol = x(3);
H2.m_sol = x(4);
M_tot_sol = 1/x(5);
disp(['solver message: ', num2str(output.message)])
disp(['No of iterations: ', num2str(output.iterations)])
disp(['The calculated minimum of Gibbs free enthalpy is: ', num2str(fval/1e9), ' GJ'])
disp(['Solution for m(CO2): ', num2str(CO2.m_sol), ' kg, solution for m(CO): ', num2str(CO.m_sol), ' kg, solution for m(H2O): ', num2str(H2O.m_sol), ' kg, solution for m(H2): ', num2str(H2.m_sol), ' kg'])

[Caroline Otto]

To my mind the problem now ist that the matrix Aeq with a determinant of -2.306956934286042e-17 is close to be singular. How can I resolve this? Do I need e.g. a QR decomposition of Aeq? If yes, how do I get the right hand side (beq)?

[Matt J]

To my mind the problem now ist that the matrix Aeq with a determinant of -2.306956934286042e-17 is close to be singular.

No, on the contrary if Aeq were a non-singular square matrix, the only feasible solution would be x=Aeq\beq, and you wouldn't need optimization at all. However, it does mean that one or more of your equality constraints may be redundant and can be removed, at least for thisnset of data.

I think it is just a matter of your initial guess, like Torsten said. It has already been demonstrated that the solution you're after is obtained once x0 is sufficiently close to it.