1
$\begingroup$

My code is based on the similarity transformation X=VZ.I simulate the model for transformed equations involving Z by replacing the state space modelAX+BU and RX+SU with transformed equations where X is replaced with VZ and WZ resp then I apply transformation again to obtain X back. I want to execute the following Euler code using solve_ivp. Although I'm able to frame the basic equations but not able to access the solution array at each time step as in Euler to obtain X back from Z. X and Z are simply the state variable matrices of order 4*1. My euler code is :

  for i in range(0,500000):
           if (w== 0 and  X[1]> vdon) or (w==1 and X[0]> 0):
                      zdot=inv(V)*A*V*Z+inv(V)*B*U
                      Z(i+1)=Z(i)+ h*zdot
                      w=1
                      X(i+1)=V*Z(i+1)
           else:
                     zdot=inv(W)*R*W*Z+inv(W)*S*U
                     Z(i+1)=Z(i)+ h*zdot
                    w=0
                    X(i+1)=W*Z(i+1)

where V and W are the eigen vector matrices of A and R respectively obtained using

e1,V=LA.eig(A)
e2,W=LA.eig(R)

In scipy's solve_ivp I define the function as

def conv(t,Z):
   if (w==0 and X[1]>vdon) or (w==1 and X[0]>0):
      zdot=inv(V)*A*V*Z+inv(V)*B*U
      w=1
   else:
       zdot=inv(W)*R*W*Z+inv(W)*S*U
       w=0
return zdot

The if condition shown here doesn't work as expected and has been shown just for code understanding. and I define the solver equations as: w=0 #intially sol= solve_ivp(conv, tspan1,X0) aa1=sol.t bb1=sol.y But I'm unable to define X=W*Z or X=V*Z at each solver time step during run time.

$\endgroup$
5
  • $\begingroup$ Could you add something about the genesis of the model. Chiefly on how the eigen-decomposition of the matrices got involved and how the unique selection of the eigen-decomposition was enforced (eigenvectors may flip signs without losing the property of being normalized eigenvectors, the order of eigenvalues might depend on the algorithm). And also add, if that is possible, a word on the "physical" meaning of Z and X. $\endgroup$ – Lutz Lehmann May 2 '20 at 10:12
  • $\begingroup$ I've edited it.Hope it brings some clarity. $\endgroup$ – hinata exc May 2 '20 at 10:19
  • $\begingroup$ You left out the phase change logic here that you had included in the previous posting of that code, that is, the changing of sw in the branches of the if condition. Could you game through and report on what should happen if X[1]>vdon is true and the second is false, that is, X[0]<0? Then you are in a state of artificial oscillation between the phases, where the oscillation frequency depends on the stepsize $h$ and not on anything internal to the model. $\endgroup$ – Lutz Lehmann May 2 '20 at 10:55
  • $\begingroup$ I added the sw changes as well.The logical operation used is or.If X[1]>vdon and sw==1 i.e. this condition is true and the other one if false, there would be no switching between the phases and the zdot=inv(V)*A*V*Z+inv(V)*B*U condition would be executed. $\endgroup$ – hinata exc May 2 '20 at 11:10
  • $\begingroup$ If sw==1, then the first clause of the first condition is false and X[1]>vdon does not get tested. Now if X[0]<0 then the second condition is false too and the else branch is entered where sw=0 is set. Now in the next step the second condition is always wrong, only X[1]>vdon gets tested, and if true then the then branch is entered and sw=1 is set. As the changes in Z and thus the computed values of X in each phase are small, this phase alternating behavior can persist for a while. This is a generalization of "sliding mode" behavior. $\endgroup$ – Lutz Lehmann May 2 '20 at 11:58
0
$\begingroup$

Independent of how that makes sense or not, you have a system with two modes of operation or phases, and the phase change condition is different for each phase, enabling something like hysteresis if the system parameters are right.

So you have the differential systems

def f0(t,Z): return A.dot(Z)+B.dot(U(t))
def f1(t,Z): return R.dot(Z)+S.dot(U(t))

and the change conditions

def g0(t,Z): return vdon(t) - V[1].dot(Z)
def g1(t,Z): return W[0].dot(Z)

The logic of your Euler code then reads as

  • if in phase 0: Integrate $Z'=f_0(t,Z)$ as long as $g_0(t,Z)>0$. Over that segment, store $X(t)=VZ(t)$. On $g_0(t,Z)=0$ switch to phase 1.
  • if in phase 1: Integrate $Z'=f_1(t,Z)$ as long as $g_1(t,Z)>0$. Over that segment, store $X(t)=WZ(t)$. On $g_1(t,Z)=0$ switch to phase 0.

This recipe has an obvious failure in that if both $g_1<0$ and $g_2<0$ then no phase is admissible, both phases are left as soon as they are entered. One could change the condition from "while $g_k>0$" to "until $g_k$ changes from positive to negative", which can be realized using events with direction

g0.terminal=True; g0.direction=-1;
g1.terminal=True; g1.direction=-1;

and then the switching mechanism is

t,z=t0,z0
sw=0
while t<t_final:
    if sw==0:
        sol = solve_ivp(f0,[t,t_final],z,events=g0)
        Z = sol.y.T
        X = V.dot(Z)
    else:
        sol = solve_ivp(f1,[t,t_final],z,events=g1)
        Z = sol.y.T
        X = W.dot(Z)
    # end if
    t,z=sol.t[-1], Z[-1]
    # collect sol.t, Z, X segments in a list to be later concatenated, 
    # see https://stackoverflow.com/a/54796425 (a discontinuous ODE entering sliding mode)
    # or directly append to large lists
#end while
$\endgroup$
3
  • $\begingroup$ Thank you for your efforts.I'll try executing the code.I didn't get the change condition properly def g0(t,Z): return vdon(t) - V[1].dot(Z) def g1(t,Z): return W[0].dot(Z) . Do you mind explaining so? $\endgroup$ – hinata exc May 2 '20 at 18:19
  • $\begingroup$ An event has to be a real-valued function. The sign of the function outside of an event is determined by the logic of the loop, not directly by the event. $\endgroup$ – Lutz Lehmann May 2 '20 at 18:43
  • $\begingroup$ ,I need to change a variable at some point t=.0005 .scicomp.stackexchange.com/questions/35063/… $\endgroup$ – hinata exc May 6 '20 at 10:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.