MATLAB: Optimization fmincon with integral. Results not optimal, what am I doing wrong

fminconintegraloptimizationOptimization Toolbox

Hi, I have the following parameterization of x,y and z

with N = 4

Then I have three forces defined as: (Clohessy-Wiltshire equations..)

and I'm trying to minimise the following cost function

which I have included in MATLAB as

% x,y,z and derivatives
x = @(t,c)  c(1) +  c(2)*t +  c(3)*t.^2 +  c(4)*t.^3 +  c(5)*t.^4;
y = @(t,c)  c(6) +  c(7)*t +  c(8)*t.^2 +  c(9)*t.^3 + c(10)*t.^4;
z = @(t,c) c(11) + c(12)*t + c(13)*t.^2 + c(14)*t.^3 + c(15)*t.^4;
xDot = @(t,c)  c(2) +  2*c(3)*t +  3*c(4)*t.^2 +  4*c(5)*t.^3;
yDot = @(t,c)  c(7) +  2*c(8)*t +  3*c(9)*t.^2 + 4*c(10)*t.^3;
zDot = @(t,c) c(12) + 2*c(13)*t + 3*c(14)*t.^2 + 4*c(15)*t.^3;
xDdot = @(t,c)  2*c(3) +  6*c(4)*t +  12*c(5)*t.^2;
yDdot = @(t,c)  2*c(8) +  6*c(9)*t + 12*c(10)*t.^2;
zDdot = @(t,c) 2*c(13) + 6*c(14)*t + 12*c(15)*t.^2;
% forces definition
uX = @(t,c) xDdot(t,c) - 3*(nOrbit^2)*x(t,c) - 2*nOrbit*yDot(t,c);
uY = @(t,c) yDdot(t,c) + 2*nOrbit*xDot(t,c);
uZ = @(t,c) zDdot(t,c) + (nOrbit^2)*z(t,c);
% cost function
fun = @(t,c) ( uX(t,c).^2 + uY(t,c).^2 + uZ(t,c).^2 ).^0.5;

The formulation of the optimization problem includes some constraints and is formulated as

f   = @(c) integral(@(t) fun(t,c),t0,tf);
x0  = ones(1,15);
A   = [];
b   = [];
Aeq = [... %Initial position and velocity constraints
    1 t0 t0^2 t0^3 t0^4 zeros(1,5) zeros(1,5); ...
    zeros(1,5) 1 t0 t0^2 t0^3 t0^4 zeros(1,5);...
    zeros(1,5) zeros(1,5) 1 t0 t0^2 t0^3 t0^4;...
    0 1 2*t0 3*t0^3 4*t0^3 zeros(1,5) zeros(1,5);...
    zeros(1,5) 0 1 2*t0 3*t0^3 4*t0^3 zeros(1,5);...
    zeros(1,5) zeros(1,5) 0 1 2*t0 3*t0^3 4*t0^3;... 
    ...%Final position and velocity constraints
    1 tf tf^2 tf^3 tf^4 zeros(1,5) zeros(1,5); ...
    zeros(1,5) 1 tf tf^2 tf^3 tf^4 zeros(1,5);...
    zeros(1,5) zeros(1,5) 1 tf tf^2 tf^3 tf^4;...
    0 1 2*tf 3*tf^3 4*tf^3 zeros(1,5) zeros(1,5);...
    zeros(1,5) 0 1 2*tf 3*tf^3 4*tf^3 zeros(1,5);...
    zeros(1,5) zeros(1,5) 0 1 2*tf 3*tf^3 4*tf^3;...
    ];
beq = [initState; finalState];
lb  = ones(1,15)*-inf;
ub  = ones(1,15)*inf;
nonlc = [];
options = optimoptions(@fmincon,'Algorithm','interior-point',...
    'MaxIterations',400,'StepTolerance',1.0e-12,...
    'MaxFunctionEvaluations',3000,'Display','notify');
s = fmincon(@(a) f(a),x0,A,b,Aeq,beq,lb,ub,nonlc,options);

I'm running this iteratively over a period of 2000 (from 500-2500)

The result I get is far from optimal and so my question is.. is there something I did wrong in the formulation of this problem? (I'm particularly unsure about the computation of the integral.) Or could it just be something with the parameters of the optimization and should I try different settings?

Resulting plot should look something like:

but instead I'm given

thank you in advance for helping out

Best Answer

The first-order optimality measure is not that small at the end. I suspect that your problem is a bit sensitive. i see in your code that you can restart the optimization using the initial parameters from the optimization previously run. That is what I was going to tell you to try, but you already knew it.

Also, your problem has several linear equality constraints. Perhaps it would be worthwhile to try the 'sqp' algorithm instead of 'interior-point' especially after you run the optimization once to get a better starting point.

Some of your parameters might be quite large or small, and perhaps the problem could use some rescaling, though I am not sure where. I mean, Opt_Testing has a mu value of about 1e14, which is fairly large. I am not sure where this parameter enters into the calculation, but I'm sure that you understand it, and might know what could be rescaled. But I am not sure that this idea has any merit.

It is possible that trying a variety of initial points could provide a better solution. It also might be worth trying a smaller value of OptimalityTolerance, perhaps 1e-10, and for the 'interior-point' algorithm try setting the ScaleProblem option to true.

Well, this is not very satisfactory advice, I am afraid, as you already seem to know very well how to use the software. Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

Related Solutions

MATLAB: Problems with fmincon, how to solve

Your problem is convex, so you can use the CVX package if you want. Also, the solution should be the global minimum; not a local minimum. To use fmincon, do the following:

Equality/nonEquality constraint:

function [c,ceq]=mycon(x)
c(1) = x(1)^2 + x(2)^2 - 4*x(1) - 3*x(2) + 5;
c(2) = -x(1)^2 - x(2)^2 + 4*x(1) + 3*x(2) - 6; % notice the fixed typo
ceq = 2*x(1)^2 + 2*x(2)^2 - 3*x(1) - 3*x(2) - 2; 
end

Construct the problem and call fmincon:

clc
clear all
Aeq=[]; 
beq=[]; 
A=[]; 
b=[]; 
% if you are not sure that vector x is within lb and ub do not restrict it
lb = []; %lb=[-3,-3]; 
ub = []; %ub=[3,3]; % 
x0=[4,4];
myObj = @(x) x(1)^2 + x(2)^2 - 5*x(1) - 6*x(2) + 15;
nonlcon = @mycon;
options = optimoptions('fmincon','Algorithm','sqp','Display','iter-detailed','MaxFunctionEvaluations',100000, 'MaxIterations',2000, 'OptimalityTolerance' ,1e-10);
[x,cost] = fmincon(myObj,x0,A,b,Aeq,beq,lb,ub,nonlcon)

Output (global minimum):


Local minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
<stopping criteria details>
x =
    1.6450    1.9007
cost =
    1.6896

MATLAB: Solve/fsolve system of equations with exp, and multiple variables.

Your d_H_tf_C(tf)==0 turns out to be a differential equation involving x(tf) . You can solve the remaining 4 equations to derive C1, C2, C3, C4, and substitute that in to d_H_tf_C(tf)==0, leaving you with a single differential equation in terms of tf and x(tf) .

When solve() sees a differential equation, it automatically invokes dsolve() on your behalf, but you could make the steps more clear by doing that solve for C1, C2, C3, C4, and trying the dsolve() on what remains.

Unfortunately what remains is too complicated for MATLAB to be able to solve.

When I pass it over to Maple, Maple generates a solution in terms of two indefinite integrals and an unresolved boundary condition. The two indefinite integrals do not appear to have closed form solutions.

EQN = -19440*exp(tf)*(-8/81*exp(3*tf)*(10*x(tf)+diff(x(tf),tf)-30)*cos(3*tf)^4+8/27*( ...
sin(3*tf)*exp(3*tf)+3/2*exp(2*tf)-3/2*exp(4*tf))*diff(x(tf),tf)*cos(3*tf)^3+((( ...
-80/9+40/27*x(tf)+4/27*diff(x(tf),tf))*exp(2*tf)-4/9*(-20+10/3*x(tf)+diff(x(tf) ...
,tf))*exp(4*tf))*sin(3*tf)+(520/81*x(tf)+16/81*diff(x(tf),tf)+200/27)*exp(3*tf) ...
+(-20/9*x(tf)+2/9*diff(x(tf),tf)-20/3)*exp(5*tf)-2/9*exp(tf)*(10*x(tf)+diff(x( ...
tf),tf)+30))*cos(3*tf)^2+2/3*diff(x(tf),tf)*((exp(tf)-22/9*exp(3*tf)+exp(5*tf)) ...
*sin(3*tf)-31/6*exp(2*tf)+31/6*exp(4*tf)-3/2*exp(6*tf)+3/2)*cos(3*tf)+((-49/27* ...
diff(x(tf),tf)-310/27*x(tf)+80/9)*exp(2*tf)+(310/27*x(tf)+25/9*diff(x(tf),tf)- ...
80/9)*exp(4*tf)+(-10/3*x(tf)-diff(x(tf),tf))*exp(6*tf)+10/3*x(tf)+1/3*diff(x(tf ...
),tf))*sin(3*tf)+(-280/27+235/81*diff(x(tf),tf)-440/81*x(tf))*exp(3*tf)+(20/3+ ...
20/9*x(tf)-29/9*diff(x(tf),tf))*exp(5*tf)+exp(7*tf)*diff(x(tf),tf)-7/9*exp(tf)* ...
(-60/7-20/7*x(tf)+diff(x(tf),tf)))/(4*exp(2*tf)*cos(3*tf)^2-22*exp(2*tf)+9*exp( ...
4*tf)+9)^2 == 0

and Maple's dsolve() says

x(tf) = (300*int(-2/27*(cos(3*tf)+1)*(cos(3*tf)-1)*(-1/3*exp(2*tf)*cos(3*tf)^2+ ...
(exp(tf)-exp(3*tf))*sin(3*tf)-7/6*exp(2*tf)+3/4*exp(4*tf)+3/4)*exp((2+3*i)*tf)/ ...
(-4/27*cos(3*tf)^4*exp(3*tf)+(4/9*sin(3*tf)*exp(3*tf)+2/3*exp(2*tf)-2/3*exp(4* ...
tf))*cos(3*tf)^3+((2/9*exp(2*tf)-2/3*exp(4*tf))*sin(3*tf)-1/3*exp(tf)+8/27*exp( ...
3*tf)+1/3*exp(5*tf))*cos(3*tf)^2+((exp(tf)-22/9*exp(3*tf)+exp(5*tf))*sin(3*tf)- ...
31/6*exp(2*tf)+31/6*exp(4*tf)-3/2*exp(6*tf)+3/2)*cos(3*tf)+(-49/18*exp(2*tf)+25 ...
/6*exp(4*tf)-3/2*exp(6*tf)+1/2)*sin(3*tf)-7/6*exp(tf)+235/54*exp(3*tf)-29/6*exp( ...
5*tf)+3/2*exp(7*tf))/((3+I)*exp((1+3*i)*tf)-4/3-i-5/3*exp(6*i*tf)),tf) + CONSTANT)* ...
exp(-10*int((2*cos(3*tf)^2*exp(tf)+(3*exp(2*tf)-3)*sin(3*tf)-2*exp(tf))/(2*cos( ...
3*tf)^2*exp(tf)+(-6*sin(3*tf)*exp(tf)+9*exp(2*tf)-9)*cos(3*tf)+(9*exp(2*tf)-3)* ...
sin(3*tf)+7*exp(tf)-9*exp(3*tf)),tf))

vpasolve() cannot deal with the system because it vpasolve() cannot deal with differential equations.

If you had one more boundary condition then you could do numeric simultations using ode45() or related (looks like it might be stiff, so perhaps ode23s). In that regard, see odeFunction()

Best Answer

Related Solutions

MATLAB: Problems with fmincon, how to solve

MATLAB: Solve/fsolve system of equations with exp, and multiple variables.

Related Question