({{, (RefLink[RiccatiSolve, paclet:ref/RiccatiSolve])[{a,b},{q,r}]
gives the matrix x that is the stabilizing solution of the continuous algebraic Riccati equation (ConjugateTranspose[a]).x+x.a-x.b.(Inverse[r]).(ConjugateTranspose[b]).x+q==0.}, {, (RefLink[RiccatiSolve, paclet:ref/RiccatiSolve])[{a,b},{q,r,p}]
solves the equation (ConjugateTranspose[a]).x+x.a-(x.b+p).(Inverse[r]).((ConjugateTranspose[b]).x+(ConjugateTranspose[p]))+q==0. }})

OpenersNotesSection(Spacer1)Details and Options

In (ConjugateTranspose[a]).x+x.a-x.b.(Inverse[r]).(ConjugateTranspose[b]).x+q==0, † denotes the conjugate transpose.

The equation (ConjugateTranspose[a]).x+x.a-x.b.(Inverse[r]).(ConjugateTranspose[b]).x+q==0 has a unique, symmetric, positive semidefinite solution x if {a,b} is stabilizable, {a,q} is detectable, r>0, and q≥0 . Consequently, all eigenvalues of the matrix a-b.(r^(-1)).b†.x are negative and the solution is stabilizing.

The solution is positive definite when {a,b} is controllable and {a,q} is observable.

(RefLink[RiccatiSolve, paclet:ref/RiccatiSolve]) supports a (RefLink[Method, paclet:ref/Method]) option, with the following possible settings:

All methods apply to approximate numeric matrices. "Eigensystem" also applies to exact and symbolic matrices.

OpenersPrimaryExamplesSection(Spacer1)Examples  (23)

OpenersExampleSection0(Spacer1)Basic Examples  (1)

Solve a continuous algebraic Riccati equation:

{a,b,q,r}={(({{-3, 2}, {1, 1}})),(({{0}, {1}})),(({{1., -1.}, {-1., 1.}})),(({{3}}))};

x=RiccatiSolve[{a,b},{q,r}]

Verify the solution:

Transpose[a].x+x.a-x.b.Inverse[r].Transpose[b].x+q//Chop

OpenersExampleSection1(Spacer1)Scope  (3)

Solve a continuous Riccati equation:

RiccatiSolve[{(({{-1, 0, 0}, {0, -4, 0}, {0, 0, -6}})),(({{3}, {1}, {-1}}))},{(({{4, 2, 0}, {2, 4, 0}, {0, 0, 1}})),{{1.}}}]//MatrixForm

Solve a Riccati equation with state-control coupling:

RiccatiSolve[{(({{0, 1, 0}, {0, -0.01, 0.3}, {0, -0.003, -10}})),(({{0, 0}, {0, -1}, {0.1, 0}}))},{(({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}})),(({{1, 0}, {0, 1}})),(({{0.1, 0}, {-0.2, 0}, {0.15, 0.2}}))}]

Solve the Riccati equation by extracting appropriate matrices from a state-space model object:

ssm=(StateSpaceModelFull[0.4245`, 0.1137`, 9.764`, 0, -0.4145`, -0.0987`, -6.32`, -3.455`, , , ]);

{q,r}={(({{84.4, 0}, {0, 6.16}})),(({{(10^(-3)), 0}, {0, (10^(-3))}}))};

RiccatiSolve[Normal[ssm][[1;;2]],{q,r}]

OpenersExampleSection2(Spacer1)Options  (7)

OpenersExampleSubsection0(Spacer1)Method  (7)

(RefLink[Automatic, paclet:ref/Automatic]) and "Eigensystem" methods can be used for exact systems:

x1=RiccatiSolve[{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},Method->Automatic]

x2=RiccatiSolve[{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},Method->"Eigensystem"]

x1==x2

As well as for inexact systems:

x1=RiccatiSolve[N@{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},Method->Automatic]

x2=RiccatiSolve[N@{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},Method->"Eigensystem"]

x1==x2

"Schur" method can be used for inexact systems:

RiccatiSolve[{N@(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},{(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))},Method->"Schur"]

"Newton" applies to inexact systems and may be more accurate than (RefLink[Automatic, paclet:ref/Automatic]):

a=SparseArray[{{1,1}->0.8,{1,5}->-0.8,{2,1}->0.8,{2,2}->1,{3,3}->-0.1,{2,5}->-0.5,{5,1}->0.25}];
b=Transpose[{{1,0,0,1,0}}];
r={{1}};
q=IdentityMatrix[5];

Compute the solution and absolute error:

x1=RiccatiSolve[{a,b},{q,r},Method->"Newton"];

Norm[(a^\[Transpose]).x1+x1.a-x1.b.Inverse[r].(b^\[Transpose]).x1+q]

Compare to the default method:

x2=RiccatiSolve[{a,b},{q,r}];

Norm[(a^\[Transpose]).x2+x2.a-x2.b.Inverse[r].(b^\[Transpose]).x2+q]

"Newton" takes suboptions "StartingMatrix", "MaxIterations", and "Tolerance":

x3=RiccatiSolve[{a,b},{q,r},Method->{"Newton","StartingMatrix"->x2,"Tolerance"->(10^(-13)),"MaxIterations"->20}];

Norm[(a^\[Transpose]).x3+x3.a-x3.b.Inverse[r].(b^\[Transpose]).x3+q]

The "Newton" method may not converge even if a stabilizing solution exists:

a=ConstantArray[0.,{5,5}];
r=b=q=IdentityMatrix[5];
RiccatiSolve[{a,b},{q,r},Method->"Newton"];

Compare with the (RefLink[Automatic, paclet:ref/Automatic]) method:

x=RiccatiSolve[{a,b},{q,r}]

"MatrixSign" is typically used as an initial approximation for the "Newton" method:

{a,b}={{{-0.1,0},{0,-0.02}},{{0.1,0},{0.001,0.01}}};
{q,r}={{{100,1000},{1000,(10^4)}},{{1.+(10^(-6)),1.},{1.,1.}}};

xinit=RiccatiSolve[{a,b},{q,r},Method->"MatrixSign"];

Use xinit to initialize the "Newton" method:

x=RiccatiSolve[{a,b},{q,r},Method->{"Newton","StartingMatrix"->xinit,"Tolerance"->(10^(-10)),"MaxIterations"->10}];

Compare errors:

Map[Norm[(a^\[Transpose]).#+#.a-#.b.Inverse[r].(b^\[Transpose]).#+q]&,{xinit,x}]

"MatrixSign" takes suboptions "MaxIterations" and "Tolerance":

x3=RiccatiSolve[{a,b},{q,r},Method->{"MatrixSign","Tolerance"->(10^(-10)),"MaxIterations"->10}];

Norm[(a^\[Transpose]).x3+x3.a-x3.b.Inverse[r].(b^\[Transpose]).x3+q]

"GeneralizedSchur" and "GeneralizedEigensystem" applies when a is singular:

a=SparseArray[{{1,1}->0.8,{1,5}->-0.8,{2,1}->0.8,{2,5}->-0.5,{2,2}->1,{3,3}->-0.1,{5,1}->0.25}];
b=Transpose[{{1,0,0,1,0}}];
r={{1}};
q=IdentityMatrix[5];

The matrix a is singular:

Det[a]

These two methods work for a singular a:

x1=RiccatiSolve[{a,b},{q,r},Method->"GeneralizedSchur"];
x2=RiccatiSolve[{a,b},{q,r},Method->"GeneralizedEigensystem"];

Verify error:

Map[Norm[(a^\[Transpose]).#+#.a-#.b.Inverse[r].(b^\[Transpose]).#+q]&,{x1,x2}]

"InverseFree" can be used when r is ill-conditioned:

{a,b}={{{-0.1,0},{0,-0.02}},{{0.1,0},{0.001,0.01}}};
{q,r}={{{10,10},{10,1000}},{{1+(10^(-9)),1-(10^(-10))},{1,1.}}};

The matrix r has a high condition number:

Norm[r] Norm[Inverse[r]]

Use the method:

x1=RiccatiSolve[{a,b},{q,r},Method->"InverseFree"]

Compare to the default method:

x2=RiccatiSolve[{a,b},{q,r}]

The absolute error is higher for the default method in this case:

{Norm[(a^\[Transpose]).x1+x1.a-x1.b.Inverse[r].(b^\[Transpose]).x1+q],Norm[(a^\[Transpose]).x2+x2.a-x2.b.Inverse[r].(b^\[Transpose]).x2+q]}

OpenersExampleSection3(Spacer1)Applications  (3)

Compute the cost associated with an optimal trajectory for a linearized inverted pendulum:

pendssm=(StateSpaceModelFull[0, 1, 0, 0, 0, ((g (m+M))/(l M)), 0, 0, 0, -(1/(l M)), 0, 0, 0, 1, 0, -((g m)/M), 0, 0, 0, (1/M), 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, , , ])/.{M-> 5.6,m-> 0.53,l-> 0.85,g-> 9.8};

{q,r}={(({{10, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 10, 0}, {0, 0, 0, 1}})),(({{1000}}))};

(Subscript[s, init])=(({{0}, {0}, {0.1}, {0}}));

The optimal cost is (Subscript[s, init]^†).x.(Subscript[s, init]):

x=RiccatiSolve[Normal[pendssm][[1;;2]],{q,r}];

(Subscript[J, opt])=(Transpose[(Subscript[s, init])].x.(Subscript[s, init]))[[1,1]]

Obtain the cost by numerical integration:

k=LQRegulatorGains[pendssm,{q,r}];

(Subscript[ss, CL])=SystemsModelStateFeedbackConnect[pendssm,k];

(Subscript[x, t])=StateResponse[{(Subscript[ss, CL]),Flatten[(Subscript[s, init])]},{0},{t,(10^5)}];

(Subscript[J, opt])=NIntegrate[(Subscript[x, t]).(q+Transpose[k].r.k).(Subscript[x, t]),{t,0,(10^5)}]

Track the "cost to go" at each instant along the optimal trajectory:

Plot[Evaluate[(Subscript[x, t]).x.(Subscript[x, t])],{t,0,50},PlotRange-> All]

Compute an optimal state-feedback gain that places all closed-loop poles to the left of -α:

{a,b,q,r}={(({{-0.2, 0.5, 0, 0, 0}, {0, -0.5, 1.6, 0, 0}, {0, 0, -14.3, 85.8, 0}, {0, 0, 0, -33.3, 100}, {0, 0, 0, 0, -10}})),(({{0}, {0}, {0}, {0}, {30}})),(({{10, 0, 0, 0, 0}, {0, 20, 0, 0, 0}, {0, 0, 6, 0, 0}, {0, 0, 0, 2, 0}, {0, 0, 0, 0, 5}})),(({{1}}))};

(Subscript[a, 1])=a +α IdentityMatrix[Length[a]]/.α-> 2;

k=Inverse[r].Transpose[b].RiccatiSolve[{(Subscript[a, 1]),b},{q,r}];

Eigenvalues[a-b.k]

The closed-loop poles without any prescribed degree of stability:

Eigenvalues[a-b.(Inverse[r].Transpose[b].RiccatiSolve[{a,b},{q,r}])]

Compute the minimum error covariance for a Kalman estimator:

{a,b,c,w,v}={(({{-1.7, 50, 260}, {0.22, -1.4, -32}, {0, 0, -12}})),(({{-272, 0.02, 0.1}, {0, -0.0035, 0.004}, {14, 0, 0}})),(({{1, 0, 0}, {0, 1, 0}})),(({{10, 0}, {0, 10}})),(({{1, 0}, {0, 1}}))};

g=b[[All,{2,3}]];

q =g.v.Transpose[g];

Tr[RiccatiSolve[{Transpose[a],Transpose[c]},{q,v}]]

OpenersExampleSection4(Spacer1)Properties & Relations  (8)

Find the optimal gains with cross-coupling matrix p:

{a,b}={(({{2.5731 +0.3269 I, 0.725 -0.3024 I}, {-0.9543-2.7184 I, -2.5731-0.3268 I}})),(({{-0.1076-0.2184I, -0.835-0.0684 I}, {-0.6379+1.198I, 0.08 +0.6334I}}))};

{q,r,p}={(({{2, 0}, {0, 4}})),(({{1, 1-I}, {1+I, 1}})),(({{4, 1}, {3, 2}}))};

x=RiccatiSolve[{a,b},{q,r,p}];

An equivalent result can be found be incorporating p into the a and q matrices:

xp=RiccatiSolve[{a-b.Inverse[r].ConjugateTranspose[p],b},{q-p.Inverse[r].ConjugateTranspose[p],r}];

xp-x//Chop

If {a,b} is controllable and {a,g} is observable, and q=(RefLink[Transpose, paclet:ref/Transpose])[g].g, then the solution to the Riccati equation is positive definite:

{a,b,g,r}={(({{0.98, 0.048, 0.0025, 0.0046, 0.008}, {0, 0.95, 0.08, 0.2, 0.533}, {0, 0, 0.23, 1.26, 6.52}, {0, 0, 0, 0.082, 1.428}, {0, 0, 0, 0, 0.367}})),(({{0.0059}, {0.509}, {10.15}, {3.97}, {1.89}})),(({{1, 0, 0, 0, 0}})),(({{1}}))};

{ControllableModelQ[StateSpaceModel[{a,b}]],ControllableModelQ[StateSpaceModel[{Transpose[a],Transpose[g]}]]}

PositiveDefiniteMatrixQ[RiccatiSolve[{a,b},{Transpose[g].g,r}]]

The eigenvalues of the Hamiltonian matrix H are pairs of the form {λ,-λ†}:

With[{a=(({{-16.826-18.09129 I, 6.1228-28.4026 I}, {19.2776 -1.14335 I, 13.826 +18.0913 I}})),b=(({{-0.8019-0.0491 I}, {0.4257 -0.4331 I}})),q=(({{10, 0}, {0, 1}})),r=(({{1}}))},H=ArrayFlatten[(({{a, -b.Inverse[r].ConjugateTranspose[b]}, {-q, -ConjugateTranspose[a]}}))]];

ListPlot[Eigenvalues[H]/.Complex[x_,y_]:>{x,y},PlotStyle-> PointSize[Medium]]

H and -H† are similar matrices:

Complement[Eigenvalues[H],Eigenvalues[-H\[ConjugateTranspose]],SameTest->(Chop[#1-#2]==0&)]

The ability to obtain a stabilizing solution depends on the Hamiltonian matrix:

{a,b,q,r}={(({{3, 0}, {0, 3}})),(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}})),(({{1, 0}, {0, 1}}))};

{vals,vecs}=Eigensystem[ArrayFlatten[(({{a, -b.Inverse[r].Transpose[b]}, {-q, -Transpose[a]}}))]]

The Hamiltonian matrix must satisfy the stability property:

Cases[vals,Complex[0,_]]

As well as the complementarity property:

stableBasis=Extract[vecs,Position[vals,_?(#<0&),{1}]];

{{(Subscript[x, 1]),(Subscript[x, 2])}}=Partition[stableBasis,{2,2}]

MatrixRank[(Subscript[x, 1])]==2

The stabilizing solution:

(Subscript[x, 2]).Inverse[(Subscript[x, 1])]//Simplify

RiccatiSolve[{a,b},{q,r}]==%

Find the eigenvalues of the feedback system with u==-(r^(-1)).(b^τ).x:

{a,b,q,r}={(({{7, 10}, {-5, -8}})),(({{(5/3)}, {-(4/3)}})),(({{0.3, 0}, {0, 0.9}})),(({{0.5}}))};

Eigenvalues[a-b.Inverse[r].Transpose[b].RiccatiSolve[{a,b},{q,r}]]

These are also the stable eigenvalues of the Hamiltonian matrix:

λ=Eigenvalues[ArrayFlatten[(({{a, -b.Inverse[r].Transpose[b]}, {-q, -Transpose[a]}}))]];

Select[λ,#<0&]

Compute optimal state feedback gains using (RefLink[RiccatiSolve, paclet:ref/RiccatiSolve]):

{a,b,q,r}={(({{0, 1}, {0, 0}})),(({{0}, {1}})),(({{1, 0}, {0, 1}})),(({{1}}))};

Inverse[r].Transpose[b].RiccatiSolve[{a,b},{q,r}]//N//Chop

Use (RefLink[LQRegulatorGains, paclet:ref/LQRegulatorGains]) to compute the same result directly:

LQRegulatorGains[StateSpaceModel[{a,b}],{q,r}]//N//Chop

Compute optimal output feedback gains using (RefLink[RiccatiSolve, paclet:ref/RiccatiSolve]):

{a,b,c,q,r}={(({{-2, 0, 1}, {0, -1, 0}, {-3, -4, -2}})),(({{0, 1}, {0, 0}, {1, 0}})),(({{1, 0, 0}, {0, 1, 0}})),(({{100, 0}, {0, 0.1}})),(({{1, 0}, {0, 1}}))};

Inverse[r].Transpose[b].RiccatiSolve[{a,b},{Transpose[c].q.c,r}]//MatrixForm

(RefLink[LQOutputRegulatorGains, paclet:ref/LQOutputRegulatorGains]) gives the same result:

LQOutputRegulatorGains[StateSpaceModel[{a,b,c}],{q,r}]//MatrixForm

Compute optimal estimator gains using (RefLink[RiccatiSolve, paclet:ref/RiccatiSolve]):

{a,b,c}={(({{-0.02, 0.005, 2.4, -32}, {-0.14, 0.44, -1.3, 30}, {0, 0.018, -1.6, 1.2}, {0, 0, 1, 0}})),(({{0.14, -0.12}, {0.36, -0.86}, {0.35, 0.009}, {0, 0}})),(({{0, 1, 0, 0}}))};

w=(10^(-2))(({{1, 0}, {0, 1}}));v=(({{(1/10000)}}));

RiccatiSolve[{Transpose[a],Transpose[c]},{b.w.Transpose[b],v}].Transpose[c].Inverse[v]

Use (RefLink[LQEstimatorGains, paclet:ref/LQEstimatorGains]) to compute this result directly:

LQEstimatorGains[StateSpaceModel[{a,b,c}],{w,v}]

OpenersExampleSection5(Spacer1)Possible Issues  (1)

If {a,b} is not stabilizable or {a,g} is not detectable, then the Riccati equation with q=(Transpose[g]).g has no stabilizing solution:

{a,b,g,r}={(({{1, 0, 0}, {1, -1, 1}, {0, 0, 2}})),(({{1}, {1}, {0}})),(({{1, 0, 0}})),(({{1}}))};

{ControllableModelQ[StateSpaceModel[{a,b}]],ControllableModelQ[StateSpaceModel[{Transpose[a],Transpose[g]}]]}

RiccatiSolve[{a,b},{Transpose[g].g,r}]

({{DocumentationFooterIconsRelatedFunction.png, ({{See Also}, {(RefLink[DiscreteRiccatiSolve, paclet:ref/DiscreteRiccatiSolve])  ▪  (RefLink[LQRegulatorGains, paclet:ref/LQRegulatorGains])  ▪  (RefLink[LQOutputRegulatorGains, paclet:ref/LQOutputRegulatorGains])  ▪  (RefLink[LQEstimatorGains, paclet:ref/LQEstimatorGains])  ▪  (RefLink[KalmanEstimator, paclet:ref/KalmanEstimator])  ▪  (RefLink[LQGRegulator, paclet:ref/LQGRegulator])}})}})

({{DocumentationFooterIconsRelatedGuide.png, ({{Related Guides}, {▪(RefLinkPlain[Control Systems, paclet:guide/ControlSystems])}, {▪(RefLinkPlain[Equation Solving, paclet:guide/EquationSolving])}})}})