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

2 ビュー (過去 30 日間)
Caroline Otto 2023 年 6 月 30 日

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
% 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

サインインしてコメントする。

### 回答 (2 件)

Torsten 2023 年 6 月 30 日

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.
##### 0 件のコメント-2 件の古いコメントを表示-2 件の古いコメントを非表示

サインインしてコメントする。

Matt J 2023 年 6 月 30 日

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;
##### 3 件のコメント1 件の古いコメントを表示1 件の古いコメントを非表示
Caroline Otto 2023 年 7 月 1 日

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 2023 年 7 月 1 日

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.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeSolver-Based Nonlinear Optimization についてさらに検索

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by