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: