1. We will be in Read Only mode (no new threads, replies, registration) for several hours as we migrate the forums to upgraded software.

[Simulink] Simulation returns to origin but it shouldn't

Discussion in 'Programmer's Corner' started by gulranek, Jun 15, 2009.

1. gulranek Thread Starter New Member

Jun 15, 2009
4
0
Hi all, this is my first post here! I posted this already on the matlab newsgroup, but no answers yet...

I'm working on a simulation in Simulink, it involves the stability of an electric power system. I'm using the Lyapunov stability equation to determine the stability of the system, dependent of two variables, current (I) and voltage (V), so we're talking about a 3-D function. The system consists of two differential equations of the first order, so I've got two integrators in Simulink (and some constants, sums etc), that's where I input my starting points. I view the V and I using the XY scope.

In my Lyapunov equation, for a certain value where V equals a constant number, let's call it n, the value of the equation is infinite, since there is a term [1/(V-n)] in the equation.

So there's basically a huge infinite wall in V=n. Basic stability theory says that if I set my initial point below V=n, there is no way my system should return to the origin (and thus to the point of stability) because there is no way it could cross over that wall (or at least not with the limited power that my system can provide).

However, I've ran dozens of simulations and no matter where I put the starting point, it always returns to the origin. I've tried reducing the fixed-step in the solver down to 10^-9, but no luck, it always crosses the barrier.

I hope I've explained myself decenlty. Does anyone have an idea as to what to do?

Thank you in advance for your suggestions!

2. steveb Senior Member

Jul 3, 2008
2,431
469
Yes, but it's too hard to go back and forth in troubleshooting Simulink problems. I've tried that before and it gets old quick.

If you are willing to post the file, I'm willing to look at it and try to figure out the issue.

3. gulranek Thread Starter New Member

Jun 15, 2009
4
0
I've sent the files to you via PM, and here's the explanation:

When z1 (voltage) equals -1, the first equation tends to infinity and the xy plot should escape to infinity, but it doesn't. My stability function consists of the two equations I've sent you, so it also goes to infinity when z1=-1.

4. steveb Senior Member

Jul 3, 2008
2,431
469
I started to look at this, but need more time. So far it looks like you've implemented the equations correctly in Simulink and I don't think there is a tolerance or solver issue. My first instinct is to say that the system is not unstable with your initial conditions you've tested so far, even if there are points or regions of instability, but I need to look at it more rigorously. It may be that initial conditions you've chosen don't go into an unstable orbit. Or maybe it is impossible to go unstable.

Clearly the equations are nonlinear so stability becomes more complex to analyze. The instability point you described seems to be only at one voltage value, so if you start at that point, you will diverge away, however, then you are no longer at that point anymore and you can eventually move into a stable region. which could put you in an orbit that attracts back to the origin.

Anyway, I need more time. The usual approach is to generate linearized state equations for the system by creating a Jacobian matrix. The matrix is dependent on the operation point (voltage and current) and the eigenvalues of the matrix (i.e. system poles) determine if that operating point is stable in a small signal (linearized) sense. Once the Jacobian matrix is identified, then we will know the unstable regions, if there are any, and whether your orbits can ever diverge away to infinity, or whether they always attract back to the origin.

One thing we can say is that you do not have the possibiliy of a choatic system because you only have two state equations. A continuous time system needs at least three states to exhibit choas. Hence I would expect the possibility of divergent solutions that go to infinity, periodic orbits and convergent solutions that converge to a point. The only way to know which are possible and from which initial conditions they occur is to map it out carefully.

Let me know if anything I've said does not seem right, but otherwise I'll go forward with this and should have an answer in about a day.

5. gulranek Thread Starter New Member

Jun 15, 2009
4
0
Yeah, everything you've said makes sense, but I did do a 3D plot of the Lyapunov function (which is a composite of the two functions I've sent you) and it has got an infinite "wall" right in V=-1, for all values of I.

So theoretically there shouldn't be a region below V=-1 where the system could return to the stable area.

Thank you again for your help!

6. steveb Senior Member

Jul 3, 2008
2,431
469
Ok, I derived the system poles and, with your conditions, the solutions are oscillatory with complex conjugate poles. This is correct since we see the spiral shape in the xy-trace. I definitely see what you are saying about the unstable points near V=-1. For the real part of the poles (Lyapunov exponent) I get the following function.

Note that once in the unstable region, the poles may not be complex conjugates all the time.

The exponent definitely becomes positive (and hence indicates instability) when near. Hence, I agree with you that this is a "wall" that you should not be able to get past. It's acutally a complex situation because if V<-V0, the rate of change is positive and the voltage is driven higher. Then, once V>-V0, the rate of change switches and the voltage is driven in the negative direction. So it's unstable, but oscillating in a sense; and perhaps mathematically the solution should stay at V=-V0, but numerically, it can't do that. This will give most numerical solutions a great deal of trouble because if the calculation is inaccurate and puts you over the wall in the stable region with positive rate of change (which happens if V>0, but can still happen if V<0 since rate of change depends on current also), you get the wrong prediction.

Looking further, I noticed that you put a saturation limit block, which gets activated in the unstable region. Once you are in saturation, your derived (i.e. linearized) equations for the exponent are no longer valid. So, with the saturation block, it may be possible to drive through the "wall", so to speak.

However, even when I disable the saturation block, a similar effect can happen, although at least the solution tries to bounce around before crossing over. My best guess is that this is just a numerical artifact that results from numerical instability.

There are several recommedations I could make.

1. Avoid the saturation blocks unless your system really is represented by that. If it is represented like that, then you must realize that the equations can't be linearized at that point where it crosses into saturation (although you can linearize the new saturated equations in a region where it's always saturated) and the Jacobian matrix representation is not valid. In this case, the numerical answers may be the only indication of the system behavior, but it's hard to be sure they are correct, because numerical instability is likely when you have positve Lyapunov exponents. This is why all weather prediction computer models give different answers, and don't always give the right answer.

2. Experiment with many solvers using the smallest timesteps and tolerances.

3. Don't use decimation on critical scopes since you may not see the real numerical behavior.

Last edited: Jun 17, 2009
7. steveb Senior Member

Jul 3, 2008
2,431
469
I was able to stabilize the solution by putting a saturation limiter block on the input of the Z1 integrator. I limited the rate of change of Z1 to plus and minus 10000 and the solution can't get over the wall. If I use a limit of plus and minus 100000, then it still does not get over the wall, but you can see numerical noise in the solution.

This may be an acceptable solution for your system because limiting rate of change on the voltage is basically a slew-rate limit and no real system is able to change voltage instantaneously.

You can also use a rate limiter on the output of the integrator, if you choose, since it has the same effect.

Last edited: Jun 17, 2009
8. gulranek Thread Starter New Member

Jun 15, 2009
4
0
Hi, sorry for the delay, I've been busy. Thanks for your help!

Here's what I managed to do using your advice: I've put the saturation limiter as you said and yes, it does stop the system from passing the wall. But it also unfortunately stops the system from returning to the origin if the starting point is in the positive part of the system (if it starts from the positive part, it should return to origin).

Does the same happen to you or not?

What type of solver and resolution did you use?

9. steveb Senior Member

Jul 3, 2008
2,431
469
I used the ODE4 (4th order Runga Kutta), but it works for many others, resolution is 1e-6 sec. Have your removed your other saturation limiter? I took that one out because I assumed that you only put it in to try and solve this issue.

Some initital conditions return to zero and some do not. For example 1,1 returns to zero, but 3,3 gets sucked into the wall. One thing about that wall at x=-V0 is that it has the ability to suck the solution into it, if the trajectory gets too close. When I started at 3,3, the solution orbits around and gets too close the the wall. In the previous post I said the following:

"It's acutally a complex situation because if V<-V0, the rate of change is positive and the voltage is driven higher. Then, once V>-V0, the rate of change switches and the voltage is driven in the negative direction."

This last sentence is relevant. Note that if V>-V0 (not over the wall yet), but is close to the wall, the rate of change becomes a very large negative number and drives the solution into the wall. It's analogous to a black-hole in space - get too close and your fall in.

Actually, this is a very interesting nonlinear system. I've never encountered a 2 state continuous system that has this feature.

As far as using the saturation limit on the state derivative. I believe this is perfectly reasonable because it equates to a slew rate limit on voltage. No doubt, your model is not including real system components and parasitics that limit slew rate. The only question is whether the 10000 V/s is consistent with your system. Note that 100000 V/s also works, but just shows numerical noise near the wall.

10. steveb Senior Member

Jul 3, 2008
2,431
469
Just to provide more information to the above post. I tried two cases and made a plot as shown.

Case 1: intial conditions Z1=Z2=1.233
Case 2: initial conditions Z1=Z2=1.234

Note that case 1 goes back to the origin and case 2 does not. This is the classic case of sensitivity to initial conditions.

File size:
15.7 KB
Views:
19