Optimization help: porting fmincon() to Julia

Specific Domains Modelling & Simulations
1

Problem

I’m having problems porting fmincon() statements that appear to require “Generic nonlinear constraints” (using nonlcon). Here are the references I have used:

Background: I just recently figured out how to port fmincon() statements successfully (I think) but I can’t quite figure out how to port “generic nonlinear constraints” like these:

options = optimset('TolX',0.001, 'TolFun',1, 'MaxIter',1000 );
options = optimset(options,'LargeScale','off');
x = fmincon(@LCObj1a,x,[],[],[],[],[],[],@LCObj1b,options,param,0.97,dbg);

With the callback functions being:

function f = LCObj1a(x,param,max_radius,dbg)
% The objective function for the initial optimization
% process used to put the roots of the denominator of the LCBP NTF inside 
% the unit circle.
f = 1;			% No actual objective
return

function [C,Ceq] = LCObj1b(x,param,max_radius,dbg)
% The constraints function for the initial optimization
% process used to put the roots of the denominator of the LCBP NTF inside 
% the unit circle.
H = LCoptparam2tf(x,param);
objective = 1;			% No actual objective
rmax = max(abs(H.p{:}));
C = rmax - max_radius;
Ceq = [];
if dbg
    fprintf(1,'x = [ ');
    fprintf(1, '%.4f ',x);
    fprintf(1,']\n');
    fprintf(1,'rmax = %f\n\n',rmax);
end
return

What I could find in Optim.jl

Under “Generic nonlinear constraints”, it appears I have to generate “generic” con_jacobian!() & con_h!() functions in addition to con_c!() (which apparently corresponds to my own LCObj1b()):

con_c!(c, x) = (c[1] = x[1]^2 + x[2]^2; c)
function con_jacobian!(J, x)
    J[1,1] = 2*x[1]
    J[1,2] = 2*x[2]
    J
end
function con_h!(h, x, λ)
    h[1,1] += λ[1]*2
    h[2,2] += λ[1]*2
end;

Then, I have to specify the following magical parameters (due to my own ignorance wrt optimization algorithms):

lx = Float64[]; ux = Float64[]
lc = [-Inf]; uc = [0.5^2]
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!,
                                     lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())

where uc should be replaced by my own radius (uc=0.97^2).

Which is great!.. except I’m still missing df. AFAICT, I have to grab the snippet from a section above:

fun(x) =  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

function fun_grad!(g, x)
g[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
g[2] = 200.0 * (x[2] - x[1]^2)
end

function fun_hess!(h, x)
h[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
h[1, 2] = -400.0 * x[1]
h[2, 1] = -400.0 * x[1]
h[2, 2] = 200.0
end;

Ok great: I have my own fun(x) (which is LCObj1a()), and I might be able to find the gradient (which is 0, isn’t it?).

Where I stopped

But without further research, I don’t really know what a Hessian is. Plus, I don’t really understand how you can converge on a solution with 0 gradient :confused: (That’s a confused face) - so I’m sort of not very confident I’m even going in the right direction.

Oh, yeah: and I’m not really sure if con_jacobian!() & con_h!() are really generic functions or if they need to be specialized for my particular application!

Help would be appreciated

Thanks in advance!

2

Hi, maybe this can help:

Basically:

  • cons is a function that,given a vector of variables,gives a vector of constraints (inplace)
  • con_Jacobian is the Jacobian of cons
  • con_hessian is a little different,it basically is an inplace adding of the Lagrange multipliers

Also, there is AD implemented for constraints, you can try using:
TwiceDifferentiableConstraints(c!, lx,ux,lc,uc,autodiff=:forward) (the documentation seems to lack behind, as the support of AD in constraints was in NLSolversBase.jl, not in Optim)

1 Like
3

@longemen3000: Your suggestions gave me enough potential leads for me to put energy into porting the code. Unfortunately, I it wasn’t sufficient for me to successfully get things working.

The code (in case you are bored or something)

→ https://github.com/ma-laforge/RSDeltaSigmaPort.jl/blob/master/src/designLCBP.jl

My problem

Part of my problem is that I don’t fully understand what the original designLCBP() function was trying to do.

  • I’m pretty certain it is trying to constrain the passband & stopband ripples to be within within a certain maximum (That’s a fairly standard thing to do).
  • I don’t fully understand the code, so I can get lost in the math (even if it is “relatively readable” Matlab code).
  • I wasn’t even aware that the pbr and sbr variables represented the passband and stopband ripples until a plot popped-up showing me what I was designing :confused: (confused).
  • I don’t really understand how the original fmincon() worked with the given inputs by reading the docs.
  • I don’t understand why I have a LCObj1a() function that has “No actual objective” (not to mention the rest of the code - which I am porting in a semi-blind fashion).

Oh yeah… and I am relatively new wrt optimization algorithms and issues related to convergence.

  • I typically use simulators. I never build them or even use optimizers in a direct fashion. SPICE-like simulators now have pretty decent presets, etc.
  • I don’t typically even need to tweak tolerances. The most I typically do is choose between TRAP & GEAR.

But thanks for having given me some leads anyhow. Now that it is 95% ported, maybe some other interested user will fix that last 5% for me :).

1 Like