How to solve implicitly constrained optimization in Matlab?

As we know, we can easily use fmincon to solve any optimization problem that contains some explicit linear and non linear, equality and inequality constraints ( such as Case 1) .

        %Case 1
         f(x) = –x1x2x3
         Sub:-x1–2x2–2x3 ≤ 0
         x1 + 2x2 + 2x3≤ 72

But, I have some additional implicit constraints , which are in a form of a nonlinear system (such as Case2) .

    %Case 2
    f(x) = –x1x2x3
    Sub: -x1–2x2–2x3 ≤ 0
          x1 + 2x2 + 2x3≤ 72
   %Here is the Additional implicit constraints in a form of nonlinear system
        X1^2+ln(u1)+sqrt(u3)=0
        X2^2+X3+ln(u2)=0
        sqrt(x1)+X3+u1=0
        X2+X1+u2^2=0
        X1^2+X2+ln(u3)=0
        X1^2+X3+ln(u1)=0
    % finally this constraints which is derived from the above system must be  %satisfied:
        0.9<real(u1,u2,u3)<1.1

As a possible but inefficient approach, I have tried to use fsolve Matlab built-in function within the constraint function file (e.g., mycon.m ) to solve the additional nonlinear system and then evaluate the added constraint 0.9<real(u1,u2,u3)<1.1 which uses new implicit variables i.e., u1,u2,u3. fimincon is a iterative method, so generally it calls mycon.m function many times (1000) to evaluate the candidate optimization vector x1,x2,x3. Thus, The nonlinear system along the added implicit constraint must be solved and evaluated as well many times. This bring a huge computational burden to this problem as fsolve is a iterative and time consuming method too.

Is this a correct approach? Is there any efficient method to include this implicit constraints 0.9<real(u1,u2,u3)<1.1 in my problem with a reasonable computational burden? Is there a way to include these implicit constraint variables (u1,u2,u3) as the optimization variables x1,x2,x3 so that fmincon evaluate them without solving the nonlinear system by fsolve? Can I consider this nonlinear system outside of 'mycon.m' so that it will be calculated once and then obtained u1,u2,u are supplied to 'mycon' for verifying the last constraint? 0.9<real(u1,u2,u3)<1.1

Note: my actual problem is huge itself and its nonlinear system contains few hundreds of equations, so using ` fsolve ` in this way, is not possible for me .

 Please tell me your suggestions.

Here I post my code for evaluating the implicit constraint using fsolve in mycon.m

%Main file
x0 = [10;10;10];    % Starting guess at the solution
[x,fval] = fmincon('myfun',x0,[],[],[],[],[],[],'mycon');

% Objective function
function f = myfun(x)
f = -x(1) * x(2) * x(3);

% Constraint function
function [c, ceq]=mycon(x)
c(1)=-x(1)-2*x(2)-2*x(3); 
c(2)=x(1) + 2*x(2 )+ 2*x(3)-72;
Y=zeros(6,1);
Y(1)=x(1);Y(2)=x(2);Y(3)=x(3);
Y(4)=1;Y(5)=1;Y(6)=1;
options = optimoptions('fsolve','Display','off');
Y = fsolve(@mysys,Y,options); % Call solver
Y=real(Y);
c(3)=Y(4)-1.1;
c(4)=Y(5)-1.1;
c(5)=Y(6)-1.1;
c(6)=-Y(4)+0.9;
c(7)=-Y(5)+0.9;
c(8)=-Y(6)+0.9;
ceq=[];

   % Nonlinear System equations.
   function eq=mysys(Y)
   x(1)=Y(1);x(2)=Y(2);x(3)=Y(3);
   u(1)=Y(4);u(2)=Y(5);u(3)=Y(6);
        eq(1)=x(1)^2+log(u(1))+sqrt(u(3));
        eq(2)=x(2)^2+x(3)+log(u(2));
        eq(3)=sqrt(x(1))+x(3)+u(1);
        eq(4)=x(2)+x(1)+u(2)^2;
        eq(5)=x(1)^2+x(2)+log(u(3));
        eq(6)=x(1)^2+x(3)+log(u(1));

[Titus Edelhofer]

Just to be sure, I understood:

• You already tried to have variables [x1, x2 x3, u1, u2, u3] instead of "only" [x1 x2 x3]?
• You are able to formulate the objective function and constraints and call fmincon?

But since your real problem has more like x1,...x100,u1,...,u100 the overall problem is to hard to solve using fmincon? Is that the challenge?

Titus

[Jamais avenir]

Dear Titus, Your first question was my problem. I didn't know that I could extend my optimization variables to [x1, x2 x3, u1, u2, u3]. Thanks for the clue.

Answer 1

[Matt J]

Is there a way to include these implicit constraint variables (u1,u2,u3) as the optimization variables x1,x2,x3 so that fmincon evaluate them without solving the nonlinear system by fsolve?

That really is what you should be doing and that's mostly all there is to say. Just rewrite your problem in terms of the longer vector of unknown design variables z=[x(:);u(:)]. This means, for example, that you'll need to add columns of zeros to your linear constraint matrices corresponding to the new variables, u. This is equivalent to saying that the u(i) don't really participate in the linear constraints.

In your mycon function, you can extract z into separate pieces x and u for your convenience. Then use the ceq output to enforce your system of nonlinear equations:

      function [c, ceq]=mycon(z)

        x=z(1:N);
        u=z(N+1:end);

            ceq(1)=x(1)^2+log(u(1))+sqrt(u(3));
            ceq(2)=x(2)^2+x(3)+log(u(2));
            ceq(3)=sqrt(x(1))+x(3)+u(1);
            ceq(4)=x(2)+x(1)+u(2)^2;
            ceq(5)=x(1)^2+x(2)+log(u(3));
            ceq(6)=x(1)^2+x(3)+log(u(1));

         c=[];

As for

     0.9<real(u1,u2,u3)<1.1

you should clarify why you are allowing u(i) to take on complex values. I'm guessing it's because fsolve isn't able to avoid complex solutions for u and that what you really want is,

      0.9<u(i)<1.1

Once you abandon fsolve, and recode everything in terms of z, you can simply enforce these bounds using the lb,ub fmincon input arguments. Note that fmincon's sqp algorithm will enforce bounds so that expressions like ln(u), sqrt(u) never give complex values. You should also be returning NaN or Inf from the objective or constraints for inputs z that don't give real and well-defined values. SQP can recover from NaNs and Infs during the search.

Finally, you should be careful with non-differentiable expressions like sqrt(u). Here, it might be okay, because you are bounding u away from u=0 where sqrt(u) is non-differentiable. Otherwise, though, you need to take differentiability issues into consideration. Options might be to rewrite

      a+b-sqrt(u)=0

as

      (a+b)^2-u=0

Or you can consider changes of variables like y=sqrt(u) so that all expressions are written smoothly in terms of y.

[Jamais avenir]

Dear Mat J Thank you for your answer. It was an illustrative example to express my problem that has such implicit constraints. However, it might not be a good example as it might not have any feasible solution at all. Also, the original last constraint can be 0.9<u(i)<1.1 as you mentioned. I put real(u1,u2,u3) because of an error in my illustrative example due to receiving complex answer by fsolve (In my earlier approach). Also, I don't want to use linear constraint matrix and want to supply them directly in nonlinear constraint file. I know it gives approximate results, but its only a example and not wrong.

Based on your suggestion I have modified my code and it works the way much faster. But, There is only a small problem in this sample problem that u(i) is out of its band [0.9,1] while I included that constraint in the 'mycon.m'. Is it because, the illustrative example that I expressed doesn't have a feasible answer, or generally its wrong to bound the variables with nonlinear inequality constraints e.g,c(6)=-u(1)+0.9;?

Here is the modified script and results:

x0 = [10;10;10;1;1;1];    % Starting guess at the solution
[x,fval] = fmincon('myfun',x0,[],[],[],[],[],[],'mycon');

   % Objective function
   function f = myfun(x)
   f = -x(1) * x(2) * x(3);

function [c, ceq]=mycon(x)
u(1)=x(4);u(2)=x(5);u(3)=x(6);
c(1)=-x(1)-2*x(2)-2*x(3); 
c(2)=x(1) + 2*x(2 )+ 2*x(3)-72;

ceq(1)=x(1)^2+log(u(1))+sqrt(u(3));
       ceq(2)=x(2)^2+x(3)+log(u(2));
        ceq(3)=sqrt(x(1))+x(3)+u(1);
        ceq(4)=x(2)+x(1)+u(2)^2;
        ceq(5)=x(1)^2+x(2)+log(u(3));
        ceq(6)=x(1)^2+x(3)+log(u(1));
        c(3)=u(1)-1.1;
        c(4)=u(2)-1.1;
        c(5)=u(3)-1.1;
        c(6)=-u(1)+0.9;
        c(7)=-u(2)+0.9;
        c(8)=-u(3)+0.9;

Results: x=
    1.7180
    3.3903
    9.0724
    0.3207
    0.0000
    0.5038

Thank you again for your helpful answer .

[Matt J]

Is it because, the illustrative example that I expressed doesn't have a feasible answer, or generally its wrong to bound the variables with nonlinear inequality constraints e.g,c(6)=-u(1)+0.9;?

It is always possible that the problem is infeasible or that you didn't choose a good enough initial guess for a feasible solution to be found.

However, lumping linear and bound constraints into the nonlcon function is also definitely not helping you. You are depriving fmincon of the knowledge that those constraints have special structure, and this makes its search more difficult, as well as more computationally demanding.

As I mentioned also, you should be using fmincon's sqp algorithm since the function and constraints here are not defined over all of R^n.

[Jamais avenir]

Thank you so much for your answer @Mat J.

Answer 2

[Torsten]

There is no other way as to use the nonlcon-function within fmincon.

Best wishes

Torsten.

[Jamais avenir]

Dear @Torsten I have edited and clarified my question. Also, Matlab script for the inefficient way that I proposed is added in the question. This code works and gives the correct result, but very slow. Please check my question again. Thank you