How to saturate a variable in ODE solver ?

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

[Torsten]

What's the error message ?

Best wishes

Torsten.

[Vincent DUFOUR]

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]

Note: This is simpler than the loop:

u = max(u, 800^2)

[Vincent DUFOUR]

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]

Sorry but doesn't change anything... I really don't get it, it's like the solver is ignoring those lines..

[Torsten]

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]

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]

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]

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]

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]

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]

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]

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]

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]

You compare the u-components to 640000, not 800.

Best wishes

Torsten.

[Vincent DUFOUR]

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

Answer 1

[Jan]

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.

[Vincent DUFOUR]

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]

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]

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]

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