Autor : Robert Muszyński KCiR PWr Pl  Wszelkie prawa zastrzeżone

Integrator Borcketta

Model

MF=MatrixForm
(eqns={ x1'[t]==u1[t],
          x2'[t]==u2[t],
          x3'[t]==x1[t]u2[t]-x2[t]u1[t]})//MF

Lagranżjan i równania Eulera-Lagrange’a

Needs["VariationalMethods`"]
sL[t_]=x1'[t]^2+x2'[t]^2 (* minimalizacja energii sterowań *)
sG[t_]=x3'[t]-x1[t]x2'[t]+x2[t]x1'[t] (* z trzeciego równania równe zero *)
sLL[t_]=sL[t]+λ[t]sG[t]
(eqEuler=EulerEquations[sLL[t],{x1[t],x2[t],x3[t]},t])//MF

Równania uproszczone

(* Równania uproszczone ze stałym λ *)
(eqEulersimp=eqEuler[[{1,2}]]/.{λ'[t]->0,λ[t]->λ}//Simplify)//MF
controls=Solve[eqns[[{1,2}]],{x1'[t],x2'[t]}][[1]]
allcontr={controls,D[controls,t]}
(eqEulersimpu=eqEulersimp/.Flatten[allcontr])//MF

(* "Rozdzielamy" zmienne *)
uprimsol=Solve[eqEulersimpu,{u1'[t],u2'[t]}][[1]]
eqsimple=D[eqEulersimpu,t]/.uprimsol//Simplify

Rozwiązanie

(usol=DSolve[eqsimple,{u1[t],u2[t]},t][[1]])//MF
(tosolve=(eqns[[{1,2}]]/.usol))//MF
(xsol= DSolve[tosolve,{x1[t],x2[t]},t][[1]])//MF

x3sol=DSolve[eqns[[3]]/.Flatten[{xsol,usol}],x3[t],t]
(xsol=Join[xsol,x3sol])//MF

(eqnfull=eqns~Join~eqsimple/.λ->2π)//MF
initcond={x1[0]==0,x2[0]==0,x3[0]==0, x1[1]==0,x2[1]==0,x3[1]==1,u1[0]==1}
finsol=NDSolve[{eqnfull,initcond},{x1,x2,x3,u1,u2},{t,0,1}];Plot[Evaluate[{x1[t],x2[t],x3[t]}/.finsol],{t,0,1},AxesLabel->{"t","x"}, PlotLegends->{"(Subscript[x, 1])","(Subscript[x, 2])","(Subscript[x, 3])"}]
ParametricPlot3D[Evaluate[{x1[t],x2[t],x3[t]}/.finsol],{t,0,1}]

Chat GPT poproszony o rozwiązanie... zgaduje jakie ono powinno być :)

(*Definiujemy równania Brocketta*)
brockettSystem={x1'[t]==u1[t],x2'[t]==u2[t],x3'[t]==x1[t] u2[t]-x2[t] u1[t]};
(*Warunki początkowe*)
initialConditions={x1[0]==0,x2[0]==0,x3[0]==0};
(*Przykładowa kontrola*)
controls={u1[t]->Cos[t],u2[t]->Sin[t]};
(*Numeryczne rozwiązanie*)
solution=NDSolve[Join[brockettSystem/.controls,initialConditions],{x1,x2,x3},{t,0,2π}];
(*Wykres trajektorii*)
Plot[Evaluate[{x1[t],x2[t],x3[t]}/.solution],{t,0,2π},PlotLegends->{"x1","x2","x3"}]

ParametricPlot3D[Evaluate[{x1[t],x2[t],x3[t]}/.solution],{t,0,2π},PlotLegends->{"x1","x2","x3"}]

A poproszony by nie zgadywał, tylko rozwiązał zadanie “naprawdę”... robi tylko błędy :( (kolejne próby po przekazaniu chatowi komunikatu o błędzie)

(*Definiujemy układ równań dla integratora Brocketta*)
brockettSystem={x1'[t]==u1[t],x2'[t]==u2[t],x3'[t]==x1[t] u2[t]-x2[t] u1[t]};

(*Warunki początkowe i końcowe*)
initialConditions={x1[0]==1,x2[0]==0,x3[0]==0};
finalConditions={x1[T]==0,x2[T]==1,x3[T]==0.5};

(*Definiujemy funkcjonał energii*)
energyCost=Integrate[u1[t]^2+u2[t]^2,{t,0,T}];

(*Rozwiązujemy problem minimalizacji energii*)
solution=NMinimize[{energyCost,Join[brockettSystem/.{x1->x1Func,x2->x2Func,x3->x3Func},initialConditions,finalConditions]},{u1,u2},Method->"Automatic"];

(*Wyświetlamy optymalne sterowanie*)
solution

Chat: Błąd w NMinimize wynika z tego, że Mathematica nie obsługuje równań różniczkowych jako bezpośrednich ograniczeń w optymalizacji. Zamiast tego należy użyć podejścia opartego na ParametricNDSolveValue, które pozwala numerycznie rozwiązać układ równań i następnie minimalizować funkcjonał kosztu.

(*Definiujemy czas symulacji*)
T=5;

(*Rozwiązujemy równania różniczkowe jako funkcje parametrów u1,u2*)
sol=ParametricNDSolveValue[{x1'[t]==u1[t],x2'[t]==u2[t],x3'[t]==x1[t] u2[t]-x2[t] u1[t],x1[0]==1,x2[0]==0,x3[0]==0},{x1,x2,x3},{t,0,T},{u1,u2}];

(*Definiujemy funkcję kosztu*)
costFunction[u1_,u2_]:=Integrate[u1[t]^2+u2[t]^2,{t,0,T}];

(*Ograniczenia na końcowe wartości stanu*)
constraints={x1/.sol[u1,u2]==0,x2/.sol[u1,u2]==1,x3/.sol[u1,u2]==0.5};

(*Minimalizacja kosztu*)
optimalControl=NMinimize[{costFunction[u1,u2],constraints},{u1,u2}];

(*Wyświetlamy wynik*)
optimalControl

Chat: Błąd wynika z tego, że ParametricNDSolveValue traktuje u1,u2u1​,u2​ jako parametry, a jednocześnie próbujemy je pozostawić jako funkcje czasu. W Mathematica sterowanie u1,u2u1​,u2​ powinno być wyrażone jako konkretne funkcje sterujące, a nie jako parametry.
Poprawne podejście:
Musimy sparametryzować sterowanie za pomocą np. funkcji interpolacyjnych lub rozwinięcia w bazie funkcji ortogonalnych (np. wielomianów). Następnie znajdziemy optymalne parametry.

(*Definiujemy czas końcowy*)T=5;

(*Definiujemy sterowanie jako wielomiany o współczynnikach {a1,a2,...,b1,b2,...}*)
u1[t_,a1_,a2_]:=a1+a2 t;
u2[t_,b1_,b2_]:=b1+b2 t;

(*Rozwiązujemy układ równań różniczkowych z danymi sterowaniami*)
sol=NDSolveValue[{x1'[t]==u1[t,a1,a2],x2'[t]==u2[t,b1,b2],x3'[t]==x1[t] u2[t,a1,a2]-x2[t] u1[t,b1,b2],x1[0]==1,x2[0]==0,x3[0]==0},{x1,x2,x3},{t,0,T}];

(*Ekstrakcja końcowych wartości stanu*)
finalState[a1_,a2_,b1_,b2_]:=Through[sol[T]];(*x1[T],x2[T],x3[T]*)(*Funkcja kosztu minimalizująca zużycie energii*)costFunction[a1_,a2_,b1_,b2_]:=Integrate[u1[t,a1,a2]^2+u2[t,b1,b2]^2,{t,0,T}];

(*Definiujemy warunki końcowe*)
constraints={First@finalState[a1,a2,b1,b2]==0,(*x1[T]==0*)Second@finalState[a1,a2,b1,b2]==1,(*x2[T]==1*)Third@finalState[a1,a2,b1,b2]==0.5 (*x3[T]==0.5*)};

(*Minimalizacja funkcji kosztu*)
optimalControl=NMinimize[{costFunction[a1,a2,b1,b2],constraints},{a1,a2,b1,b2}];

(*Wyświetlamy wynik*)
optimalControl

Chat: Błąd NDSolveValue::ndnum oznacza, że Mathematica napotkała wartości nienumeryczne dla pochodnej, ponieważ a1, a2, b1, b2a1​, a2​, b1​, b2​ są symboliczne, a NDSolveValue wymaga konkretnych liczb.🔧 
Poprawne rozwiązanie:
Musimy wymusić, aby NDSolveValue otrzymywało numeryczne wartości sterowania.Zrobimy to poprzez SetDelayed ( := ) oraz Evaluate. (co oczywiście nic nie pomogło)

(*Definiujemy czas końcowy*)T=5;

(*Układ równań różniczkowych*)
eqns={x1'[t]==u1[t],x2'[t]==u2[t],x3'[t]==x1[t] u2[t]-x2[t] u1[t],λ1'[t]==-λ3[t] u2[t],λ2'[t]==λ3[t] u1[t],λ3'[t]==0,(*λ3 to stała*)(*Warunki brzegowe*)x1[0]==1,x2[0]==0,x3[0]==0,x1[T]==0,x2[T]==1,x3[T]==0.5,(*Optymalne sterowanie*)u1[t]==-1/2 (λ1[t]-λ3[t] x2[t]),u2[t]==-1/2 (λ2[t]+λ3[t] x1[t])};

(*Rozwiązujemy układ*)
sol=NDSolveValue[eqns,{x1,x2,x3,λ1,λ2,λ3,u1,u2},{t,0,T}];

(*Wyciągamy rozwiązania*)
{x1Sol,x2Sol,x3Sol,u1Sol,u2Sol}=sol[[{1,2,3,7,8}]];

(*Rysujemy wykresy trajektorii*)
Plot[{x1Sol[t],x2Sol[t],x3Sol[t]},{t,0,T},PlotLabels->{"x1","x2","x3"}]

(*Wykresy sterowania*)
Plot[{u1Sol[t],u2Sol[t]},{t,0,T},PlotLabels->{"u1","u2"}]

I w końcu poproszony by rozwiązał zadanie stosując równania Eulera-Lagrange’a i metodę mnożników Lagrange’a... dalej robi tylko błędy :(

(*Definiujemy czas końcowy*)
T=5;

(*Układ równań różniczkowych*)
eqns={x1'[t]==u1[t],x2'[t]==u2[t],x3'[t]==x1[t] u2[t]-x2[t] u1[t],λ1'[t]==-λ3 u2[t],λ2'[t]==λ3 u1[t],λ3'[t]==0,(*λ3 to stała*)(*Warunki brzegowe*)x1[0]==1,x2[0]==0,x3[0]==0,x1[T]==0,x2[T]==1,x3[T]==0.5,(*Optymalne sterowanie*)u1[t]==-1/2 (λ1[t]-λ3 x2[t]),u2[t]==-1/2 (λ2[t]+λ3 x1[t])};

(*Rozwiązujemy układ*)
sol=NDSolveValue[eqns,{x1,x2,x3,λ1,λ2,λ3,u1,u2},{t,0,T}];

(*Wyciągamy rozwiązania*)
{x1Sol,x2Sol,x3Sol,u1Sol,u2Sol}=sol[[{1,2,3,7,8}]];

(*Rysujemy wykresy trajektorii*)
Plot[{x1Sol[t],x2Sol[t],x3Sol[t]},{t,0,T},PlotLabels->{"x1","x2","x3"}]

(*Wykresy sterowania*)
Plot[{u1Sol[t],u2Sol[t]},{t,0,T},PlotLabels->{"u1","u2"}]