MATLAB: How to solve implicitly constrained optimization in Matlab

fminconfsolvenonlinearoptimizationOptimization Toolbox

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));

Best Answer

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.
Related Question