How to saturate a variable in ODE solver ?

2017 3 15
1 个回答
更新时间：2017 3 16
7 次查看（30 天）

0 个投票

Hi everyone,
I am currently using ode15s to solve a ODE system and I would need to saturate a variable directly in the dynamic, is it possible ?
To make it clear, I am resolving the following system:
if true
p-dot-dot=(1/M)*wRb*F1*u+f
w_dot = inv(J)*F2*u+g
end
wRb, J , F1 and F2 are matrices, f and g other terms and I would need to prevent my "u" variable from being more than a threshold.
I already tried to do this:
for i=1:length(u)
if u(i)>800^2
u(i)=800^2;
end
end
just after the computation of u and before computing p_dot_dot and w_dot, unfortunately it doesn't seem to work...
Thanks in advance

16 个评论
显示 14更早的评论 隐藏 14更早的评论

Torsten
Torsten 2017-3-15
What's the error message ?
Best wishes
Torsten.
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
There is no error, when I compare the result with my loop inside the ode solver and without, my 'u' is the same.
I even tried to put this loop in a other function that is called to change the value of 'u' but even this doesn't change anything..
Here are my file, don't bother with all the commented code you can just run the mainODE. What I need is my "wi" variable in the main to be saturated at 800 (the 'u' variable is wi^2).
Thanks
Jan
Jan 2017-3-15
Note: This is simpler than the loop:
u = max(u, 800^2)
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Max will replace all the value less than 800 of u by 800 I want the inverse, but I think using min it should work.
Thanks for the idea Jan I would have never thought about it
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Sorry but doesn't change anything... I really don't get it, it's like the solver is ignoring those lines..
Torsten
Torsten 2017-3-15
编辑：Torsten 2017-3-15
If u is a numerical input array to a function, it cannot be modified.
For the local use within the function, a copy of u is made, and this copy can be modified.
http://blogs.mathworks.com/loren/2006/05/10/memory-management-for-functions-and-variables/
Best wishes
Torsten.
Jan
Jan 2017-3-15
Yes, of course, min() would be correct. Sorry. This was thought only to let the code run (or crash) faster and more elegant. The problem is somewhere else. By the way: You did not explain yet, what the problem is. I do not know, what "u" is also. The solver does not ignore anything, it does not have the power to do so. Please edit the question and add more details.
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
U is not a numerical input, it is computed directly by the solver then used to compute the solution of the equation and I just want to saturate it that's all I don't get why I can't..
Thanks a lot
Torsten
Torsten 2017-3-15
As I told you: uin can not be modified.
Use
function u = compar(uin,threshold)
u=uin;
for i=1:6
if uin(i,1)>threshold
u(i,1)=threshold;
end
end
return
Best wishes
Torsten.
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Oh sorry Torsten I misunderstood I thought you were talking about the u to be an input of the solver.
I tried your way and still have the same u not being modified..
Torsten
Torsten 2017-3-15
u will not be modified, but the output u_sat in the calling program will contain the saturated version of u.
Best wishes
Torsten.
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Yeah sure, sorry I wasn't clear, I meant that I put my exit dX(8:13,1)=sum(abs(u_sat)) but still I received the same value for my u-sat than for my u..
I suppose I have to head to the event function direction to see what I can do.
Thanks again
Torsten
Torsten 2017-3-15
编辑：Torsten 2017-3-15
Maybe all u-components are smaller than 800^2 :-)
Or you must call "compare" like
threshold=640000;
u_sat=compare(u,threshold);
Or the following bug-report is relevant for your case:
https://de.mathworks.com/matlabcentral/answers/102610-why-do-global-variables-not-update-in-the-workspace-browser-or-the-array-editor-when-i-modify-them-f
Best wishes
Torsten.
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
No they are I am trying to do it because they are bigger than 800 and that's not realistic.
Thanks for your time trying to help
Torsten
Torsten 2017-3-15
You compare the u-components to 640000, not 800.
Best wishes
Torsten.
Vincent DUFOUR
Vincent DUFOUR 2017-3-16
All right guys to conclude the topic I actually realized that I was using the wrong data... I was extracting the 'u' through dX=u but it logically made the sum over the time which doesn't mean anything for me here...
The saturation actually works !
Thanks all of you for your time and help
Cheers

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Jan
Jan 2017-3-15

1 个投票

Matlab's ODE integrators are designed to handle functions with a ontinuous derivative. Functions like max, min, if, abs etc. makes the output non-smooth and in consequence they should not appear inside the function to be integrated. The resulting discontinuities confuse the step size control. The results beein used to calculate one step can be coming from intervals with different definitions for "u". If you are lucky you get an error message, but without luck Matlab squeezes the stepsize until the rounding error dominates the final value and the trajectory is more or less random.
See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047.
Please do not post "doesn't seem to work..." without providing any details. My answer is based on guessing what going wrong.

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Oh okay I didn't know about the min and max function but you're the one who proposed it at first ^^
"Doesn't seem to work" just to mean it doesn't change anything, but as I can't find any answer I suppose it's not possible using solver..
Thanks for all the answers
Jan
Jan 2017-3-15
It can be implemented using event function: Whenever u reachs a limit, the integration is stopped and restarted using a different parameter for the function to be integrated.
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Ok thanks I'll see in this direction, btw I uploaded my files to be run if you want to see.
To make it quick, it's the ODE system of a six-rotor drone, u=wi^2 and wi are the rotation speed of the propeller which I want to be constrained (<800rad/s), the p_dot_dot is the acceleration in position and w_dot in rotation.
In this case I want my drone to follow a sinus trajectory on pitch and roll, while all the other (x,y,z and yaw) remain at 0.
Thanks a lot for all your answers
Vincent DUFOUR
Vincent DUFOUR 2017-3-15
Other question for Jan, is it possible to set an event on my variable 'u' which is not my state x ?
I saw that to use the Event option of odeset the function will have the form [val,isterminal,direction] = myEventFcn( t,x) but I need an event on my 'u' not x.. don't know if it's possible

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by