MATLAB Answers

How to get fmincon to work with an ode45?

5 views (last 30 days)
JV
JV on 17 Nov 2020
Commented: Star Strider on 17 Nov 2020
Hello fellow members of the community
I've written a code for the purpose of finding the specific initial values of an ODE of thrown particle. I wanted it be based on ode45 and fmincon function.
So much to say, it is not working and spits out an error
I'll be very thankful for any help
Cheers
Jan
PS
Here I've copied full code, that should work on any computer.
fun = @Throw;
x0 = [10 6];
nonlcon = @Ceilling;
[v0opt, fval] = fmincon(fun,x0,[],[],[],[],[],[],nonlcon);
function dpos = Throw(t,pos)
mu = 0.0094;
g = 10;
dpos = zeros(4,1);
dpos(1) = pos(3);
dpos(2) = pos(4);
dpos(3) = -mu * pos(3)* sqrt(pos(3)^2 + pos(4)^2);
dpos(4) = -g - mu * pos(4) * sqrt(pos(3)^2 + pos(4)^2);
end
function j = FindParam(v0)
targetPosZ = 5;
targetPosY = 2;
Opt = odeset('Events', @StopIntegration);
[t,pos] = ode45(@Throw,[0, 6],[0, 1, v0(1), v0(2)],Opt);
plane = plot(pos(:,1),pos(:,2),[targetPosZ, targetPosZ],[0 3.5],targetPosZ,targetPosY,'o');
j = (targetPosY - pos(end,2))^2 % obj. fun. = (dist. to target)^2
end
function [c,ceq] = Ceilling(v0)
ceq = [];
targetPosZ = 5;
targetPosY = 2;
Opt = odeset('Events', @StopIntegration);
[t,pos] = ode45(@Throw,[0, 10],[0, 1, v0(1), v0(2)],Opt);
h = v0(2)^2/2*g;
v = sqrt(v0(2)^2 + v0(1)^2);
c(1) = max(h - 5); % height of the ceilling
c(2) = min(y(2) + 0); % projectile can't touch the floor
c(3) = max(v - 20); % max velocity of 20
end
function [value, isterminal, direction] = StopIntegration(t, y)
targetPosZ = 5; % integration stops when projectile hits the wall
value = (y(1) > targetPosZ);
isterminal = 1;
direction = 0;
end

  0 Comments

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 17 Nov 2020

  2 Comments

JV
JV on 17 Nov 2020
I found it - I had used wrong function handle. Still, I will for sure look at that documentation, so thanks!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by