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: