MATLAB Answers

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

53 views (last 30 days)
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

  2 Comments

Alan Weiss
Alan Weiss on 15 Dec 2020
That's a pretty impressive problem formulation. Can you please give us the exit flag, function value, and output structure when you run
[s,fval,exitflag,output] = fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlc,options);
I could not run the code because I did not have the t0 input for the A matrix.
Alan Weiss
MATLAB mathematical toolbox documentation
Nico de Korte
Nico de Korte on 16 Dec 2020
The output I get from the (first) run is this:
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>
fval:
1.6318
exitflag:
1
output:
iterations: 9
funcCount: 174
constrviolation: 4.4107e-07
stepsize: 0.1631
algorithm: 'interior-point'
firstorderopt: 0.0819
cgiterations: 27
message: '↵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>↵↵Optimization completed: The relative first-order optimality measure, 7.533789e-13,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 1.128680e-20, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
I've included my test files as .zip so you can also have a look for yourself if you want to.

Sign in to comment.

Accepted Answer

Alan Weiss
Alan Weiss on 16 Dec 2020
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

  3 Comments

Nico de Korte
Nico de Korte on 16 Dec 2020
Many thanks for your answer. I will definitely try some of the options and working with the sqp. I was already afraid of this sensitivity, but I guess it will have to be dealt with manually then.
David Goodmanson
David Goodmanson on 16 Dec 2020
Hello NIco,
It appears that you are specifying both the initial and final positions. I take it that are you specifying the time from initial to final position and integrating over that to find the cost, is that correct? Or is the integration time an arbitrary time somewhere along the path?
Could you provide some typical values for initial and final time, position, and velocity, and nOrbit? In other words, everything it takes to run the problem, including anything else I might have forgot?
Nico de Korte
Nico de Korte on 16 Dec 2020
Hi David,
If you see my answer to Alan's earlier question you will find, in an attached .zip file, my code which is plug-and-play and where all values that you are asking for are provided.
Indeed I am specifying both initial and final positions as well as velocities. This gives a total of 4*3 = 12 constraints on x,y,z. Since N = 4, x,y and z are parameterized by 5 coefficients and thus a total of (3*5 =) 15 coefficients is estimated, meaning (15-12=) 3 free coefficients for the optimization.
The loop runs from 0 to 2500 with update intervals of 500 seconds, with the first 500 seconds 'free motion'. Thus optimization from 500-> 2500, 1000->2500 ... etc ... 2000->2500.
Please have a look at the .zip file for any coefficients you might need. If you have further questions, let me know. Thanks in advance for your time.

Sign in to comment.

More Answers (0)

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!