Lagrange multipliers in constrained nonlinear optimization

I want to find the minimum of the following constrained nonlinear model using Lagrange multipliers.

I have written code from internet and modified it according to my model as follows :

clc
clear all


syms x1 x2 x3 x4 x5 x6 x7 x8 miu1 miu2 psi1 psi2 psi3 psi4 s0 s1 s2
s0 = x1^2 + x2^2 + x3^2 + x6^2 - 4*x1 - 6*x2 - 8*x3 - 62*x6 + 626;
s1 = x2^2 + x3^2 + x5^2 + x6^2 + x8^2 - 6*x2 - 8*x3 - 38*x5 - 62*x6 - 14*x8 + 719;
s2 = x2^2 + x3^2 + x4^2 + x7^2 + x8^2 - 6*x2 - 8*x3 - 14*x4 - 20*x7 - 14*x8 + 151;
A = x1^2 - 4*x1*x3 - x2*x6 + 1.57*s0;
B = x5*x6 + x8*x3 - x2*x3 + 0.71*s1 - 240 == 0;
C = 2*x7*x3 + x8*x3 - x2*x4 + 1.12*s2 - 70 == 0;
D = x1*x8 - x5 == 0;
E = x3*x4 - x6 == 0;
F = x2^2 - x7 == 0;
G = x1*x2 - x8 == 0;


L = A + miu1 * lhs(B) + miu2 * lhs(C) + psi1 * lhs(D) + psi2 * lhs(E) + psi3 * lhs(F) + psi4 * lhs(G); % Langrange function

dL_dx1 = diff(L,x1) == 0; % derivative of L with respect to x1
dL_dx2 = diff(L,x2) == 0; % derivative of L with respect to x2
dL_dx3 = diff(L,x3) == 0; % derivative of L with respect to x3
dL_dx4 = diff(L,x4) == 0; % derivative of L with respect to x4
dL_dx5 = diff(L,x5) == 0; % derivative of L with respect to x5
dL_dx6 = diff(L,x6) == 0; % derivative of L with respect to x6
dL_dx7 = diff(L,x7) == 0; % derivative of L with respect to x7
dL_dx8 = diff(L,x8) == 0; % derivative of L with respect to x8

dL_dmiu1 = diff(L,miu1) == 0; % derivative of L with respect to miu1
dL_dmiu2 = diff(L,miu2) == 0; % derivative of L with respect to miu2
dL_dpsi1 = diff(L,psi1) == 0; % derivative of L with respect to psi1
dL_dpsi2 = diff(L,psi2) == 0; % derivative of L with respect to psi2
dL_dpsi3 = diff(L,psi3) == 0; % derivative of L with respect to psi3
dL_dpsi4 = diff(L,psi4) == 0; % derivative of L with respect to psi4

system = [dL_dx1; dL_dx2; dL_dx3; dL_dx4; dL_dx5; dL_dx6; dL_dx7; dL_dx8; dL_dmiu1 ; dL_dmiu2; dL_dpsi1; dL_dpsi2; dL_dpsi3; dL_dpsi4]; % build the system of equations
[x1_val, x2_val, x3_val, x4_val, x5_val, x6_val, x7_val, x8_val, miu1_val, miu2_val, psi1_val, psi2_val, psi3_val, psi4_val] = solve(system, [x1 x2 x3 x4 x5 x6 x7 x8 miu1 miu2 psi1 psi2 psi3 psi4], 'Real', true); % solve the system of equations and display the results 
results_numeric = double([x1_val, x2_val, x3_val, x4_val, x5_val, x6_val, x7_val, x8_val, miu1_val, miu2_val, psi1_val, psi2_val, psi3_val, psi4_val]); % show results in a vector of data type double

when the code is run, the result is not output.

which part of the code needs to be fixed? Or is the whole code wrong?

[John D'Errico]

Why do you assume there must be an algebraic solution to the problem? You have what is effectively a moderately high order polynomial problem. (It is that. Consider that if you use those equaity constraints, to substitute one variable for another, thus reducing the number of unknowns, as you do so, you increase the order of the exponents on the variables? All of those products and squares of variables combine, into what is effectively a high degree polynomial problem, at least if you could do all of the algebra.) Do you understand that general polynomials of degree higher than 4 will have no algebraic solution?

Solve is a symbolic solver. It tries to find an analytical solution to the problem, in the form of radicals, etc. So solve will generally fail here. The presience of inequality constraints means there may be infinitely many solutions, and that makes it more problematic yet for solve, since solve looks for a unique solution.

The point is, you are using the wrong tool to solve a problem of this complexity. Symbolic tools are generally not good at problems like this. Even if a solution could be found, the solution is usually far more time sonsuming than you would tolerate. Instead, you need to use a numerical solver. Of course, if you do, then you really don't need to use Lagrange multipliers, since a tool like fmincon does that work for you.

But, could you use Lagrange multipliers to solve the problem, perhaps using fsolve? That should be entirely possible. Is this homework? Is there a valid reason why you need to do this the hard way?

Answer 1

[John D'Errico]

Let me give a simple example, showing why solve is a poor tool for such a problem.

minimize x^2 + y^2 + z^2

subject to x+x*y+z == 1

Not even an inequality. Solving it via Lagrange multiplers should be simple, right? It should be a fairly well-posed problem. The minimum if I have no constraint must lie at (0,0,0). And the point at the origin does not satisfy the equality constraint. So there must be some well defined minimum.

syms x y z lambda real
L = x^2 + y^2 + z^2 + lambda*(x + x*y + z - 1)

L = 

dL = gradient(L)

dL = 

Note there is no need to differentiate every term. JUST USE GRADIENT!

The solution now requires a solution of that set of 4 equations in 4 unknowns. Gosh, this problem is simple!

sol = solve(dL == 0,[x,y,z,lambda],'returnconditions',true)

sol = struct with fields:
             x: - u^3/8 - u^2/4 + u/2 + 1
             y: u^4/16 + u^3/8 - u^2/4 - u/2
             z: -u/2
        lambda: u
    parameters: u
    conditions: u^5 + 2*u^4 - 8*u^3 - 16*u^2 + 32*u + 32 == 0 & in(u, 'real')

So there is an unknown parameter u introduced. And all we need to know is the algebraic solution to that 5th degree polynomial in u. We can find a numerical solution. But solve must fail. Why? BECAUSE ALGEBRAIC SOLUTIONS TO THAT 5TH DEGREE POLYNOMIAL CANNOT BE FOUND. This goes back to long before we were born, to the days of Abel-Ruffini.

solve(sol.conditions)

ans = 

You could use a numerical solver.

sol = vpasolve(dL == 0,[x,y,z,lambda])

sol = struct with fields:
         x: 0.49170929568498349342332969358149
         y: 0.20126887244695420513765473761056
         z: 0.40932492880081387410505947868792
    lambda: -0.81864985760162774821011895737585

Is that the same solution that fmincon would find?

X = optimvar("X",[1,3]);
prob = optimproblem;
prob.Objective = sum(X.^2);
prob.Constraints.C1 = X(1) + X(1)*X(2) + X(3) == 1;
X0.X = [2 2 2];
solve(prob,X0)

Solving problem using fmincon.

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.

ans = struct with fields:
    X: [0.4917 0.2013 0.4093]

And of course fmincon found the same solution.

But you need to understand that even for an almost trivial problem to formulate as I posed above, solve failed to yield a solution, because the problem reduces to a polynomial of degree 5.

Now go back and look at the problem you have posed? What do you expect will happen?

Could we have posed the same problem to fsolve? Surely, yes.

fsolveobj = @(xyzL) [xyzL(1)+xyzL(3)+xyzL(1)*xyzL(2)-1; ...
    2*xyzL(1)+xyzL(4)*(xyzL(2)+1); ...
    2*xyzL(2)+xyzL(4)*xyzL(1); ...
    xyzL(4)+2*xyzL(3)];

[xyzL,fval,exitflag] = fsolve(fsolveobj,[2 2 2 2])

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

xyzL = 1×4

0.4917 0.2013 0.4093 -0.8186

fval = 4×1

1.0e-10 * 0.1204 0.1004 0.0598 0

exitflag = 1

And again, fsolve has no problem solving it, just as I expect, fsolve could in theory solve your problem too. Certainly, if a solution exists and the problem is well-posed, fmincon should succeed.