Sequential Convex Programming

[Johannes 1]

Hi Everyone,

I have a problem applying "sequential convex programming" to a jerk constrained path tracking optimizer.

The code i have written is based on the paper, that i attached. My guess is, that i have set my constraints poorly within the outer iteration loop, dealing with the linearized problem.

I have already introduced some slack sdp variables, as it is recomended by Prof. Löfberg for debugging, without any success. Also commenting constraints led to no result.

I hope someone can help me with this problem and thanks in advance.

Greets,

Johannes

[Johan Löfberg]

If path variable -s- unequal distributed, check all numerical differentiation!!

Undefined function or variable 'derive1'.

[Johannes 1]

Oh sorry, i forgot to attach these functions!! 

[Johan Löfberg]

it's pretty trivially infeasible

optimize(constraint(9))

K>> optimize(constraint(9))

CPXPARAM_Read_APIEncoding                        "*"

CPXPARAM_MIP_Display                             1

Infeasibility or unboundedness row 'c1':  ax <= -1.12382e+042.

Presolve time = 0.00 sec. (0.00 ticks)

ans = 

  struct with fields:

    yalmiptime: 1.0067

    solvertime: 0.0683

          info: 'Infeasible problem (CPLEX-IBM)'

       problem: 1

Judging from what cplex tells us, you should look at your data is it currently involves pure numerical crap.

[Johannes 1]

Ok, i have found some mistakes in the constraint formulation and now it runs through. The only crux i still struggle with is setting bvar(1) and bvar(K+1) to zero without running into numerical difficulties...

Nevertheless, thanks a lot Johan!!

[Johan Löfberg]

You probably end up with ill-formed SOCPS etc. 

For instance

constraint = [constraint, cone([2*cvar(1,k); bvar(1,k)-1], bvar(1,k)+1)];

will say cvar^2<=0 when bvar is 0, which perhaps might lead to poor geometry of the SOCP. Maybe you should just add the equality cvar==0 in those cases

[Johannes 1]

That was the point, now it runs just perfect!!!

Thanks a lot for the hint Johan!!

[Johan Löfberg]

BTW, you appear to do a bunch modelling to encode sum(cpower(bvar),-0.5)

...which I would change to sum(cpower(bvar(2:end-1)) as you force bvar(1) and bvar(end) to be 0, so that discrete approximation will be infinite

[Johan Löfberg]

However, as always when I see all these iterative convex schemes, the important question is how well it performs compared to simply solving the original NLP using an NLP solver. In the end, any high-level iterative convex approximation scheme to solve an NLP, is just another high-level implementation of an NLP solver, which can be extremely costly.

[Johan Löfberg]

should be sum(cpower(bvar,-0.5))

[Johannes 1]

Ok, this sounds logical!! So the next step I will take is to built up my problem by using  sum(cpower(bvar,-0.5))... If i will get this running the next step for me will be the use of a NLP solver instead of the iterative approach and check the performance of the different methods!

Again thanks a lot for your support Johan, it's a really great help!!

[Johan Löfberg]

Yes, that would be interesting (I have an old publication myself on socp-based path optimization, so it would be interesting to see  straight NLP performance)

[Johannes 1]

Ok great, i'll let you know if i have some results!!

[Johannes 1]

Ok, so the sequential approach runs with 'CPLEX', but not with 'SEDUMI' or 'SDPT3', which I can't understand.

Your proposal of skipping the sequential approach and formulate a non-convex problem was not solveable, maybe I formulated the problem wrong. 

So as a mather of fact, I will stick to the sequential approach and will try to figure out, why 'SEDUMI' of 'SDPT3' can't solve the problem.

[Johan Löfberg]

if cplex works then sedumi and sdpt3 should work, unless you have integer variables (or a nonconvex model where integers are introduced to model the nonconvex use of operators)

I would be interested in the code for the nonlinear approach, to see if anything can be done

[Johannes 1]

Hi Johan! I retried the nonconvex approach and now it runs through... It took FMINCON 500 Iterations (max Iterations was reached i guess) to solve the problem and the result i got was not really "optimal"... By introducing a slack - variable for the jerk constraint, the solver doesn't converge at all.

I attached the code snipped!

So far, i think the sequential convex approach is the best choice for the problem! It just takes a few seconds to find the optimum for a 4 dimensional problem and the result i get is really optimal in terms of exploiting the consttraints. Nevertheless, i still need to have a closer look, on why SEDUMI and SDPT3 don't work here.

[Johan Löfberg]

so what does the code look like for the version which sdpt3/sedumi isn't applicable?

[Johan Löfberg]

A problem is that you are still using the cone operators, which really isn't relevant once you go for a general nonlinear solver (and I wouldn't be surprised if there are bugs in the gradient callbacks for inequalities using cone in a nonlinear solver, as that is a very little tested feature)

Hence, simply using 

constraint = [constraint, cvar(1,k)^2 <= bvar(1,k)];

and

constraint = [constraint, 1 <= (cvar(k+1)+cvar(k))*dvar(k), dvar(k)>=0];

is a more natural model for a nonlinear solver

fmincon has some issues (interior-point solver appears to find a reasonable solution rather quickly but does not terminate, some setting I guess. sqp solver has problems with an initial point where the derivative of the square root is infinite, can probably be fixed)

However, ipopt solves the problem in 0.3 seconds

[Johannes 1]

Sorry for the delay, I had to make up a code snipped for the sequential approach first. And the whole sequential approach (building and solving) takes 16.5 seconds which is of orders higher than the nonconvex one by using IPOPT, as you suggested. So I guess that's the method of choice!!

Yet I first need to get IPOPT on my machine, which seems to be a bit of a heck, but nothing to resign over i guess.... 

[Johan Löfberg]

On my machine it takes 25 seconds of which only 2.5 seconds is spent in the solver. Hence, efficient use of the optimizer framework would probably take this down to 3s or so, as it looks like your model is easily parameterized in a few changing elements. Still smoked by ipopt though 

[Johannes 1]

Hi Johan, i still struggle with IPOPT!

When I include YALMIP and OPTI into my path, so that i can call IPOPT as a solver, i get the following error:

Attempt to execute SCRIPT ipopt as a function:

C:\Users\herbj\Documents\MATLAB\OPTIMIZATION\MODELLING\OPTI-master\Solvers\ipopt.m

Error in callipopt (line 122)

[xout,info] = ipopt(model.x0,funcs,options);

Error in solvesdp (line 350)

    eval(['output = ' solver.call '(interfacedata);']);

Error in optimize (line 31)

[varargout{1:nargout}] = solvesdp(varargin{:});

Error in minimal_nonconvex>path_tracking_optimizer_jerk_nonlin (line 207)

sol_diagnostics = optimize (constraint,objective,options);

Error in minimal_nonconvex (line 52)

[ time_opt, qs_opt, avar, bvar, jerk_opt_var, accel_opt_var, vel_opt_var ] = path_tracking_optimizer_jerk_nonlin...

Since OPTI includes a precompiled MATLAB interface, my opinion was, that i can call IPOPT just like any other solver... 

[Johan Löfberg]

This looks like what happens when there is a missing binary. So unless you've forgot to include some subdiredtory, you don't have the correct binaries

On the last version of opti, I had to run the installation scipt to install all binaries

[Johan Löfberg]

>> which ipopt

C:\work\solvers\OPTI-Master\Solvers\ipopt.mexw64

[Johannes 1]

Hi Johan,

thanks a lot for the tip! Yet i had some problems including the Intel MKL. 

After reading a bit in the OPTI-Forum, i purchased another way to get the latest .MEX files. I ran the opti_Install.m file, which guides you to the latest .MEX files on github. Now i can include IPOPT with no effort...

[Johannes 1]

It seems like the reason for the sequential approach to fail with sedumi are the numerically calculated derivatives of the sdp variable bvar... 

Is it delicate to include derivatives of sdp variables in optimization problems?

[Johan Löfberg]

not sure what you mean as conic solvers never make callbacks for derivatives

[Johannes 1]

For the jerk constraint, it is necessary to include the derivative of the sdp variable bvar with respect to the path increment ds. 

I calculated them numerically and set them to the constraints:

    %------------------------------------------------------------------------------------------------------%

    % first derivative of bvar: twopoint formula

    for k=1:K+1

        if k==K+1

            constraint = [constraint, [(0-bvar(k))/ds(k-1) == dbvar(k)]:'Derivative last']; % endpoint

        elseif k==1

            constraint = [constraint, [(bvar(k+1)-0)/ds(k) == dbvar(k)]:'Derivative first']; % startpoint

        else

            constraint = [constraint, [(bvar(k+1)-bvar(k))/ds(k) == dbvar(k)]:'Derivative mid']; % startpoint

        end

    end

    %------------------------------------------------------------------------------------------------------%

    

    %------------------------------------------------------------------------------------------------------%

    % second derivative of bvar: threepoint fromula

    for k=1:K+1

        if k==1

            constraint = [constraint, ...

                [( 2*bvar(k)/(ds(k)*(ds(k)+ds(k+1))) -2*bvar(k+1)/(ds(k)*ds(k+1)) +2*bvar(k+2)/(ds(k+1)*(ds(k)+ds(k+1))) )...

                == ddbvar(k)]: 'DDerivative first']; % startpoint

        elseif k>=2 && k<=K

            constraint = [constraint, ...

                [( 2*bvar(k-1)/(ds(k-1)*(ds(k-1)+ds(k))) -2*bvar(k)/(ds(k-1)*ds(k))  +2*bvar(k+1)/(ds(k)*(ds(k-1)+ds(k))) )...

                == ddbvar(k)]: 'DDerivative mid']; % midpoin

        elseif k==K+1

            constraint = [constraint, ...

                [( 2*bvar(k-2)/(ds(k-2)*(ds(k-2)+ds(k-1))) -2*bvar(k-1)/(ds(k-2)*ds(k-1)) +2*bvar(k)/(ds(k-1)*(ds(k-1)+ds(k-2))) )...

                == ddbvar(k)]: 'DDerivative last']; % endpoint

        end

    end

    %------------------------------------------------------------------------------------------------------%

    

I also tried to introduce slackvariables for the derivative constraints, which wasn't successfull...

[Johan Löfberg]

so you simply get an ill-conditioned problem for some reason

start could be to write T*df = f1-f2 instead of df = (f1-f2)/T, and similiar equivalent things