Optimize a group of parameters to minimize the delay which is constrained by a second-oder ODE, PLEASE HELP!

The problem is detailed in the attachment.

In short, I have a second-order ODE. No closed form solution can be found. Unfortunately, the objective function I need to optimize IS depending on the solution of the ODE. In addition, there are also some other equations/inequalities which set more constraints. The good side is there equations/inequalities are analytically given. How should I opmtimize this problem? Which Matlab function should I use? Please help. Thanks a lot!

[RahulTandon]

MAY I PLEASE SOLVE THE ODE FOR YOU? MIGHT HELP IN PART IF NOT FULLY! [e-mail] .

Answer 1

[Alan Weiss]

Without reading your problem, perhaps you can find some help in this example and in this information.

In short, I would try to use fmincon or fminunc to do the optimization, but would carefully try to choose an appropriate value of 'FinDiffRelStep' or 'DiffMinChange'.

And if you need to pass extra parameters to the optimization, see Passing Extra Parameters.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Chuang Qian]

Thanks! But could you explain the first equation in this example ? Why is the ODE written in that way?

[Alan Weiss]

I am not sure what you are asking. The first equation is a 2-component vector equation for the second derivative (force) of a body, subject to the force of gravity and to a quadratic force in the negative direction of the current velocity.

For the second equation, the original second-order ODE is turned into a 4-component first-order ODE system using the usual procedure to reduce higher-order ODEs to first order systems.

Alan Weiss

MATLAB mathematical toolbox documentation

[Chuang Qian]

What does ||x(t)||x(t) mean?

To my understanding,the motion equation should be like the following: (assume the horizontal and vertical coordinates are x1 and x2 respectively, )

d(dx1(t)/dt)/dt=-0.02v1(v1^2+v2^2)^0.5 (horizontal)

and d(dx2(t)/dt)/dt=-0.02v2(v1^2+v2^2)^0.5-9.81 (vertical),

where v1 and v2 are the velocity in horizontal and vertical direction. So I think the left-hand side of this equation should be

||dx(t)/dt||dx(t)/dt, rather than ||x(t)||x(t).

Do I make any mistakes here?

[Alan Weiss]

You are quite correct, I made a mistake in the original documentation equation. Indeed, the equation should be

d^2x/dt^2 = -0.02||v(t)||v(t) - (0,9.81),

where v(t) = dx/dt.

You will notice that the MATLAB code in fact implements that correct equation, as x_3(t) and x_4(t) are the velocity components.

I will correct the documentation. Thanks for pointing out my mistake.

Alan Weiss

MATLAB mathematical toolbox documentation

[Chuang Qian]

First of all, I should thank you for the suggestion. Things are close to work, but not exactly there yet.

I tried both 'fmincon' and 'patternsearch'. Both seemed running, but neither honored the constraint 'ceq'. I wished the error of 'ceq' to be within 1e-18,and I set 'TolCon' to 1e-19. But both 'patternsearch' and 'fmincon' gave ceq(x_opt) to about 1e-16, which was way too high. Any suggestions to improve this? Thanks a lot!

[Alan Weiss]

Your problem most likely needs to be rescaled. In floating point arithmetic, the relative error of a calculation can be no better than about 2e-16, and is often worse. So if you are talking about a relative error of 1e-18, then you are out of luck using floating point arithmetic.

But if you need your absolute error to be smaller than 1e-18, then you need to scale your problem so that this absolute error corresponds to a relative error of 1e-12 or so. In other words, multiply your ceq function by 1e6 or some such thing. You might need to rerun your optimization after it finishes.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 2

[RahulTandon]

    %%declarations , you can change here
syms Y em b k E A t gzero gd k Fa Vdd clear;
em = 1.23e-16;
b = 0;                      % for simplicity
E = 8.85e-12;
A = 1e-12;
gd = 5e-9;
Fa = 1.8e-9;
Vdd = 1;                   % for ease!!
gzero = 1;                 %%???!!!
k = 0.363;                  %%ha ha !
clc;

%% symbolic evaluation Eqn1 = 'em*D2Y+b*DY+k*Y == (E*A/(2*(gzero-Y)))*(Vdd+(Y*(k*gd-Fa)*(gzero-gd)^2)^(1/2))^2'; Eqn2 = 'Y(0) == 0'; Eqn3 = 'DY(0) == 0'; % Y = dsolve(Eqn1,Eqn2,Eqn3); % did not solve on m PC, i shud try parallel % pooling etc. %% next, get the function handle M = matlabFunction(em*diff(Y,t,2)+b*diff(Y,t)+k*Y - (E*A/(2*(gzero-Y)))*(Vdd+(Y*(k*gd-Fa)*(gzero-gd)^2)^(1/2))^2,'vars',{'t' 'Y'});

%% third, , grafical/numerical evaluation clc; [T,Y] = ode45(M,[-10 10],[0 0]); clf; plot(T,Y);