Message Boards Message Boards

0
|
9859 Views
|
1 Reply
|
0 Total Likes
View groups...
Share
Share this post:

How to solve algebraic looping in Mathematica ?

Posted 9 years ago

Hi guys so here is the code

(*Parameters*)

tend = 100;
Mu = 1;
Rho = 2.62*10^-12 Exp[-(500000 - r[t] )/64500];
Rtet = Rearth + 500000 - r[t] ;
ugc = 6.67*10^11; 
massEarth = 5.9721986*10^24; 
v = Sqrt[(ugc*massEarth)/Rtet]
c1 = 2.67;
a1 = 8;
m1 = 1000;

c2 = 2.5;
a2 = 4;
m2 = 500;

B1 = c1*a1/m1;
B2 = c2*a2/m2;

hh = 64500;
dt = 0.005;

k1 = 1;
k2 = 1;

(*Disturbance*)

Q1 = v^2*r[t]*Rho*Cos[Theta]Cos[Phi]*((Mu*(B1-B2)/2)-   
     (dt*hh)*(1-Cos^2[Phi[t]]Sin^2[Theta[t]])^0.5/Cos[Theta[t]]Cos[Phi[t]]);

Q2 = v^2*r[t]*Rho*Sin[Theta]Sin[Phi]*((Mu*(B2-B1)/2)-   
     (dt*hh)*(1-Cos^2[Phi[t]]Sin^2[Theta[t]])^0.5/Cos[Theta[t]]Cos[Phi[t]]);

(*Flip Flop Control Law*)

 u1 = k1*Theta'[t]-(Q2+u2)*(Theta'[t]+1)/Phi'[t];

 u2 = k2*Phi'[t]-(Q1+u1)*Phi'[t]/(Theta'[t]+1);


(*Anti Disturbance Control Law*)

sigma = (u1+Q1/(2*Theta'[t]+1))+(u2+Q2/(2*Phi'[t]));


(*System Dynamics*)

e1 = Theta''[t]+(Theta'[t]+1)*(2*sigma-       
2*Phi'[t]*Tan[Phi[t]])+3*Sin[Theta[t]]*Cos[Theta[t]]-Q1;

e2 = Phi''[t]+2*sigma*Phi'[t]+
((Theta'[t]+1)^2+3*Cos^2[Theta[t]])*Sin[Phi[t]]*Cos[Phi[t]]-Q2;

e3 = r'[t]-sigma*r[t];


(*Initial conditions*)

system = NDSolve[{e1 == 0, e2 == 0, e3 == 0, Theta[0] == 0.1, 
Theta'[0] == 0, Phi[0] == 0.1, Phi'[0] == 0.001, r[0] == 1}, {Theta,Phi,r},
{t,0,tend}, MaxSteps -> Infinity];


(*Plot*)

Plot[Evaluate[Theta[t] /.system], {t,0,tend}, Frame -> True, ImageSize ->  
{400, 300}, PlotStyle -> Red, LabelStyle -> Directive[12], FrameLabel ->  
{"time(s)", "Theta (rad)"}, PlotRange -> All]

Plot[Evaluate[Phi[t] /.system], {t,0,tend}, Frame -> True, ImageSize ->  
{400, 300}, PlotStyle -> Blue, LabelStyle -> Directive[12], FrameLabel ->  
{"time(s)", "Phi (rad)"}, PlotRange -> All]

So the problem here lies in the (Flip Flop Control Law ) section. Since the equation in u1 is dependent on u2 and vice-versa, therefore this results in an algebraic looping and one method I know of in order to overcome this problem is to introduce a unit delay. In SIMULINK, the "Unit Delay" block is available, but in Mathematica, how do I program a unit delay for this system so I can simulate and plot the results successfully for Theta and Phi ? Thank you so much for your effort and time. I truly appreciate it and attached is a Mathematica program file to this problem for your reference.

Best regards, Aaron Aw

Attachments:
POSTED BY: Aaron Aw
  1. The first step to debugging this and many other problems is to create as simple of an example of the problem as possible. Can you write an example of a simple differential equation which shows this problem?

  2. Create Functions. They're the key to good programming for people new to Mathematica. Both Q1, Q2, u1 and u2 are functions that depend on t. So you'd write them as:

    u1[t_] := k1Theta'[t] - (Q2[t] + u2[t])(Theta'[t] + 1)/Phi'[t];

    u2[t_] := k2Phi'[t] - (Q1[t] + u1[t])Phi'[t]/(Theta'[t] + 1);

Start with this. Turn everything that is a function into a function.

POSTED BY: Sean Clarke
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract