Converting nonlcon Constraints into YALMIP Constaints in Automated Non-Black Box Way

[Mark L. Stone]

Suppose I have nonlinear constraints specified in a nonlcon constraints function in FMINCON form.  There may or may nor be gradients provided as optional 3rd and 4th arguments, which I suppose YALMIP would be free to ignore in favor of its own gradient computation.

My paradigm is to use this within my own MATLAB function which, in order o solve subproblems, either calls a solver directly (i.e., without using YALMIP)  or via YALMIP.  The function handle of the nonlinear constraints function is an input to my MATLAB function, and is passed though to any solver called directly, and that works fine.  I would like to add a capability to allow solvers to be called via YALMIP, making use of  the nonlcon function handle.  Assume that all the code in the nonlcon function to create a c (inequality with rhs = 0) or ceq (equality to 0) would be valid YALMIP code if entered in the normal way in YALMIP.

Example:

x = sdpvar(4,1)
Constraints = [my bound and linear constraints];
Constraints = [Constraints, <nonlinear constraints in nonlcon function handle>]

% Here is an example nonlcon function
function [c, ceq, GC, GCeq] = nonlcon_example(x,p)
c = [x(1)^2+x(4)^2-2; x(2)^3+2*x(3)^5-3-p];
ceq =  [];
if nargout==4
GC = [2*x(1) 0 0 2*x(4);0 3*x(2)^2 10*x(3)^4 0]';
GCeq = [];
end


So I want to add c <= 0 to my bound and linear constraints.  How do I do this in a "good" way, so that for instance, YALMIP can compute the constraint gradient and provide it to the solver?  I.e., how do I create  <nonlinear constraints in nonlcon function handle>

Thanks.

[Johan Löfberg]

Not sure what you are asking here Mark. Did you accidentally copy text from KiDo Tom or something, because this sounds like something you know very well already?

[Mark L. Stone]

No, this is my own question which I've had on the back burner for a while.  I have a special purpose optimizer still under development for which I currently have FMINCON, KNITRO, and BARON as solvers available for subproblems having nonlinear constraints.  I just pass though the nonlcon function handle (with adjustment to the bounds if calling BARON in order to match its nonlinear constraint conventions).  I call those solvers directly without using YALMIP.  I would like to add other solvers, such as PENLAB, in order to support LMIs and BMIs in addition to nonlcon style nonlinear constraints. I would like to be able to use any existing nonlcon function which is already directly suitable for use in FMINCON (or KNITRO, which is the same).

Incorporating bound (lb and ub) and linear constraints (A, Aeq, b, beq) into my YALMIP problem formulation is simple. I am trying to figure out how nonlinear constraints from fmincon style nonlcon function handle can be inserted into my YALMIP constraints. I know exactly how to do this if I sit at the keyboard and enter it; I am trying to determine how this would be automated for use inside my optimizer. I would then add LMIs and BMIs, if any, which I know how to do, into the YALMIP constraints. I declare my YALMIP variables, such as  x = sdpvar(n) .  Then what?

Thanks.

[Johan Löfberg]

Not sure what you are asking. Perhaps sdpfun is what you are looking for?

http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.Sdpfun

[Mark L. Stone]

No, sdpfun is not what I want. 

I want to get a nonlinear constraint, consisting entirely of functions and operations which are fully supported in YALMIP, into my YALMIP formulation. This would be extremely straightforward if the constraint was manually entered by me into a YALMIP "session".  The difficulty for me, I am guessing/hoping not for you, is because my nonlinear constraint(s) are specified inside of a FMINCON style nonlinear constraints function nonlcon, for which I have the function handle available.  Inside my MATLAB function (top level optimization function), I wish to set up a YALMIP problem formulation which includes bound, linear, and perhaps LMI/BMI constraints. I know how to do that and call a solver such as PENLAB or maybe another solver, using optimize and sdpsettings. I don't know how to get the nonlinear constraints out of the function specified by the function handle, and added to my other YALMIP constraints.  So I wish to "import" the nonlinear constraints in the function handle into a YALMIP constraint, in a "good" way, such that YALMIP can provide gradient of my constraints, rather than treating them like a black box.

Going back to the example in my first post, and letting p be a parameter passed in the function handle, I wish to "extract"  [x(1)^2+x(4)^2-2; x(2)^3+2*x(3)^5-3-p]  <= 0 and add it to my constraints which I will use in the optimize command.  But I am not there to copy and paste it in; this has to occur in an automated way inside my top level optimization function.

Thanks.

[Johan Löfberg]

Why can you not simply apply the function?

f = @(x)sin(x)
 sdpvar x
 f(x)
 Nonlinear scalar (real, models 'sin', 1 variable)

[Mark L. Stone]

I wish to use existing nonlinear constraints functions, developed for use in FMINCON or KNITRO, in YALMIP.  Some of the nonlinear constraints functions could be considerably more complicated than the example.  When my top level function is called, the nonlinear constraints function handle is an input. 
So in the example, @(X)nonlcon_example(X,3) woudl be one of the input arguments (parameter p = 3 in this invocation).  This will be read in as nonlcon in my top level optimization function.  I wish to have the constraints in nonlcon extracted and placed into (added to) YALMIP constraints, with which I can call optimize with my choice of supported solver.  My top level optimization function provides the objective function and bound, linear and LMI/BMI constraints, if any, and I know how to enter all of these into the YALMIP optimize formulation, but I do not how to enter the nonlinear constraints in nonlcon.

[Johan Löfberg]

I don't see how one would be able to extract the model from a black-box function handle

[Mark L. Stone]

Oh well, now you can understand why i didn't konw how to do it.  But I was hoping you had some kinf of variable extract, string, eval kid of combination or something which would do it. I was hoping you would know how to see inside the black box.

[Mark L. Stone]

O.k., to change the "rules" a bit.  Suppose that the file handle for ny nonlinear constraints function is made available inside my program, as well as the handle (with parameters) used to invoke it.  So in this example, @(X)nonlcon_example(X,3) as well as the file handle for the m file containing nonlcon_example are available.  Then could the constraints be extracted?  You can assume that my function has already done a feval on this function handle, and therefore knows which of c (1st output argument which is for <= 0 constraints) and ceq (2nd output argument which is for ceq == 0 constraints). is nonempty, and therefore wishes to extract the constraints out of whichever or both of c and ceq are non-empty.

[Johan Löfberg]

No idea

[Mark L. Stone]

I'll look into mtree. Thanks.  I need to find some documentation.

[Johan Löfberg]

mtree perhaps? Never used it and have no idea how it works

parsetree = mtree('array2table.m', '-file')

parsetree = 

  mtree (complete: 180 nodes)

>> 

>> help mtree

 mtree  Create and manipulate M parse trees

    This is an experimental program whose behavior and interface is likely

    to change in the future.

[Xiaolong Jin]

Hi Mark, I am also being haunted by the problem... I just have no idea about converting nonlcon constraints function in FMINCON form into YALMIP constaints, since the function handle of the nonlinear constraints function is an input to my MATLAB function.
So could you kindly give me some suggestions about this matter?
Thank you very much!

[Johan Löfberg]

Perhaps you should take a step back, and describe what kind of model/constraint you want to model using YALMIP, putting your previous fmincon implementation aside

[Xiaolong Jin]

Thank you for your kind reply

I define the "RCS" as a binary variable to represent the operation status of a electric distribution network, with 1 for closed and 0 for open.

RCS=binvar(1, 37); %%%The IEEE 33 feeder has 37 branches

The constraint for "RCS" is to keep the electric distribution network a radial structure, so I have to call a radial check function as follows:

F1=[@(RCS)radial_check(RCS)==1]; %%%if the constraint for electric distribution network is satisfied, the function return value of "radial_check" would be "1"

If the constraint for electric distribution network is satisfied, the function return value of "radial_check" would be "1"

But it does not work by calling the "radial_check" to express the constraint "F1"...

I know the objective function handle can be conducted by using the "dspfun" command in YALMIP just like:

sdpvar x 
y = sdpfun(10*x,'sin')
optimize([-pi <= x <= pi],y+abs(x))

So, is there a command can be called to model the constraint by function handle just like the objective?

Could you kindly give me some suggestions?

[Johan Löfberg]

You don't want to use such a nasty non-differentiable impulse function.

What does the radial check look like?

For instance, if it checks if x'*x <=1, means you would write a simple convex well-defined constraint as a non-differentiable discontinuous black-box object, which of course would be utter madness (just an example, but almost all cases would be much better do directly, instead of converting it to a logical expression through a black box)

[Xiaolong Jin]

Thank you very much for your kind suggestions

I will have a try to model the constraint directly anyway...

Thank you again for the amazing YALMIP, which really helps me a lot in my research

[Tom]

Hi,

I was wondering if there are any new solution approaches concerning this problem (considering it's been almost three years since the last post).

Is there any approach which can be considered optimal for implementing such constraints in YALMIP, or should I turn to developing my own Matlab RO code. Are there any experiences with this isssue...

Any help would be appreciated.

Best regards

[Johan Löfberg]

not sure I understand what you want to be able to do (perhaps as I don't understand what RO stands for in this context, and the fact that the old thread appears to involve multiple questions about different things)

[Tom]

The problem that I am trying to solve is basically very similar to the original one posted by Mark.

I have an optimization algorithm already coded in Matlab which I use for power system optimization. Optimization model includes also nonlinear constraint functions. In order to calculate these constraints I have to precalculate some other functions (using the provided iteration values of x), which are not directly part of the constraint set. For that I use the nonlcon Matlab function (link), in which I define all nonlinear constraints.

My idea was to implement some uncertainties in the existing optimization model using the YALMIP robust optimization (that's where the RO comes from, sorry for being unclear).

So, the question is are there some ways to implement the defined nonlcon function in the YALMIP model (as well as the existing objective function)?

Thanks for the help.

[Johan Löfberg]

OK, no support for black-box  operators (beyond sdpfun)

[Tom]

Ok. Thank you for your quick response.

[Mark L. Stone]

matbar, which is the interface of BARON to MATLAB, and was programmed by Jonathan Currie, the OPTI Ttoolbox developer,, can read in a handle for a nonlcon style nonlinear constraints function, and get the nonlinear constraints from the nonlcon file into the internal format needed by BARON.  So I would think  think that a similar approach, whatever that is, could work for YALMIP.

help baron
 baron Solve a NLP/MINLP using the global MINLP solver baron
   
    Calling Form:
        [x,fval,exitflag,info,allsol] = baron(fun,A,rl,ru,lb,ub,nlcon,cl,cu,xtype,x0,opts)

...

 x = baron(fun,...,ub,nlcon,cl,cu) solves subject to the nonlinear
    constraints specified by nlcon, cl and cu of the form
    cl <= nlcon(x) <= cu. nlcon is a column-wise vector function and must
    be a function specified using baron rules. To specify a one
    sided constraint use infinity on the empty side, and cl = cu for an
    equality constraint.

...

  Matlab functions supplied to this method are processed and passed to
    baron as a text model and therefore must meet certain requirements. An
    example function is shown below:
 
        obj = @(x) 3*x(1)^2 - 2*x(2)^3 + 3*x(1)*x(2);
 
        OR
 
        function fx = objective(x)
           fx = 3*x(1)^2 - 2*x(2)^3 + 3*x(1)*x(2);
        end
 
    The function must be:
        - A function of one input argument (the decision variable vector)
        - Contain scalar, vector, or matrix operations only (max 2D supported)
        - Limited to a subset of functions (exp, log, log10, dot, sqrt, norm, sum, prod)
        - Contain no conditional statements or relational operators with respect to the decision variables
 
    For nonlinear constraints the function must be a column-wise vector:
 
        con = @(x) [-(x(1)*x(2));
                    3*x(1) + 2*x(2)];
 
        OR
 
        function cx = nlconstraints(x)
           cx(1) = -(x(1)*x(2));
           cx(2) = 3*x(1) + 2*x(2);
        end

[Ankit Singh]

Hello Mark,

actually, I am also kind of stuck with the same issue as yours where I have to convert the nonlcon constraints into YALMIP. Did you have any breakthrough with this?

Kind regards,

Ankit

[Johan Löfberg]

What is it you actually want to do?

This is no problems, as long as those expressions are built using YALMIP-compatible expressions

x = sdpvar(4,1)
Constraints =[];

[c, ceq] = nonlcon_example(x,p)

Constraints = [Constraints,c,ceq];

[Mark L. Stone]

You forgot the `<= 0 ` and =`= 0`  I.e., should be

Constraints = [Constraints,c <= 0,ceq == 0];

and that works.

This solution is so simple and straightforward, it's hard for me to imagine it eluded me and (apparently, unless you misunderstood the question), you. Was there a change to YALMP which makes this possible now, but not in the YALMIP of 5 years go?

[Johan Löfberg]

Right.

Yes, looks like we were confused back then. I had some memory of the issue being getting the derivatives into yalmip, but reading it now that does not seem to have been the issue

[Ankit Singh]

Dear Prof. Löfberg,

I have created a reproducible non-linear optimization example which uses 'fmincon' as a solver, but it seems there is a compatibility issue (by the way I am solving this on OCTAVE):

Here is the code:

x = sdpvar(3,1);

obj = -x(1)*x(2) - x(2)*x(3);

C = [x(1)^2 - x(2)^2 + x(3)^2 <= 2];

C = C + [x(1)^2 + x(2)^2 + x(3)^2 <= 10];

assign(x, [1; 1; 0]);
ops = sdpsettings('usex0', 1, 'debug', 'on', 'solver', 'fmincon');
sol = optimize(C, obj);

The output of 'sol' is:

sol =

  scalar structure containing the fields:

    yalmipversion = 20200116
    yalmiptime =  NaN
    solvertime =  NaN
    info = Unknown problem in solver (Turn on 'debug' in sdpsettings) (Invalid call to fmincon.  Correct usage is:

 -- Function File: fmincon (OBJF, X0)
 -- Function File: fmincon (OBJF, X0, A, B)
 -- Function File: fmincon (OBJF, X0, A, B, AEQ, BEQ)
 -- Function File: fmincon (OBJF, X0, A, B, AEQ, BEQ, LB, UB)
 -- Function File: fmincon (OBJF, X0, A, B, AEQ, BEQ, LB, UB, NONLCON)
 -- Function File: fmincon (OBJF, X0, A, B, AEQ, BEQ, LB, UB, NONLCON,
          OPTIONS)
 -- Function File: fmincon (PROBLEM)
 -- Function File: [X, FVAL, CVG, OUTP] = fmincon (...))
    problem =  9

What do you think? Am I missing something?

Kind regards,

Ankit

[Johan Löfberg]

You should have debug set to 1 to begin with

[Ankit Singh]

But, even then the result is still the same.

[Johan Löfberg]

but it looks like fmincon in octave doesn't support a trailing param argument, making it incompatible with fmincon

[Johan Löfberg]

if octave supports lambda functions try this change in callfmincon

if model.linearconstraints
    callback_con = [];
else
    callback_con = @(x)fmincon_con_liftlayer(x,model);
end

solvertime = tic;
f = @(x)fmincon_fun_liftlayer(x,model);
[xout,fmin,flag,output,lambda] = fmincon(f,model.x0,model.A,model.b,model.Aeq,model.beq,model.lb,model.ub,g,model.options.fmincon);

[Ankit Singh]

[xout,fmin,flag,output,lambda] = fmincon(f,model.x0,model.A,model.b,model.Aeq,model.beq,model.lb,model.ub,g,model.options.fmincon);

I assume, in this code 'g' means 'callback_con', right?

[Johan Löfberg]

yes typo

if model.linearconstraints

    g = [];

else

    g = @(x)fmincon_con_liftlayer(x,model);

[Ankit Singh]

Hi,

successful update:

I made the changes as you suggested, however since 'lambda' is not supported in Octave so I removed it from the return list and you made a small mistake with 'g' in the fmincon arguments, so I set it to callback_con. So now it works wonderfully.

Thanks Professor.

[Johan Löfberg]

Ok, I've taken those fixes and pushed to the develop branch. Thanks for checking.

[Ankit Singh]

Hello everyone,

this morning, I installed the new version of YALMIP to try out new updates. However, I am getting a new kind of error for my nonlinear optimization problem which is:

'newEqualities' undefined near line 33 column 17

But I can't figure out where is this error originating.

Could someone help?

Thanks

[Johan Löfberg]

Well I guess you found the first bug. You ill have to supply a reproducible example.

[Ankit Singh]

Although in this reproducible example, I am not getting the above mentioned error, but I am still getting another error:

Example:

x = sdpvar(3,1);

obj = -x(1)*x(2) - x(2)*x(3);

C= [];

C = [x(1)^2 - x(2)^2 + x(3)^2 <= 2];
C = C + [x(1)^2 + x(2)^2 + x(3)^2 <= 10];
assign(x, [1; 1; 0]);

ops = sdpsettings('usex0', 1, 'debug', 1, 'solver', 'fmincon');
sol = optimize(C, obj, ops);
value(x)

Error:

error: cat: dimension mismatch

[Johan Löfberg]

I cannot reproduce that.

Version of MATLAB? OS?

[Ankit Singh]

Its in Octave: version
ver
----------------------------------------------------------------------
GNU Octave Version: 5.2.0 (hg id: eb46a9f47164)
GNU Octave License: GNU General Public License
Operating System: MINGW32_NT-6.2 Windows 6.2  x86_64

----------------------------------------------------------------------

OS: Windows 10

[Johan Löfberg]

Ah. Well then it is completely unrelated to the previous error you reported. I found that, and that can only be in the solution of a mixed-integer problem

Line 8 in presolveTrivialSDP should be moved outside the loop to line 6.

[Johan Löfberg]

I've fied some other isues and pushed, but that particular error does not seem to be a YALMIP issue, but something more intricate inside fmincon which happens for some iterates. You will have to debug that 

error: cat: dimension mismatch

error: called from

    __constraints_interface__>@<anonymous> at line 107 column 2

    __lm_feasible__ at line 156 column 13

    fmincon at line 417 column 24

    callfmincon at line 95 column 29

    solvesdp at line 368 column 5

    optimize at line 31 column 24

debug> backend (f.objf, pin, hook)

debug> dct = df_cstr (p, ac_idx, setfield (dfdp_hook, "f", v_cstr))

error: cat: dimension mismatch

error: called from

    __constraints_interface__>@<anonymous> at line 107 column 2

    __lm_feasible__ at line 156 column 13

    fmincon at line 417 column 24

    callfmincon at line 95 column 29

    solvesdp at line 368 column 5

    optimize at line 31 column 24

debug>

[Johan Löfberg]

BTW, I noticed either my Octave or general Octave is completely broken. sdpsettings takes forever to complete here (it should take 0.1s or so)

tic;ops = sdpsettings('usex0', 1, 'debug', 1, 'solver', 'fmincon');toc

Elapsed time is 31.6697 seconds.
>>

what do you have on your installation

[Ankit Singh]

 tic;ops = sdpsettings('usex0', 1, 'debug', 1, 'solver', 'fmincon');toc

Elapsed time is 2.22305 seconds.
>>

It works fine for me it seems. I just have the GNU Octave installed as usual. No idea why it is so.

[Johan Löfberg]

OK better (but still way too slow compared to MATLAB)