Wolfram Language & System Documentation Center

DiscreteRiccatiSolve

gives the matrix that is the stabilizing solution of the discrete algebraic Riccati equation .

solves .

Details and Options

In , denotes the conjugate transpose.
The equation has a unique, symmetric, positive semidefinite solution only if is stabilizable, is detectable, , and . Consequently, all the eigenvalues of the matrix lie inside the unit circle, and the solution is stabilizing.
The solution is positive definite when is controllable and is observable.
DiscreteRiccatiSolve supports a Method option. The following settings can be specified:

	Automatic	automatically determined method
	"Eigensystem"	based on eigen decomposition
	"GeneralizedEigensystem"	based on generalized eigen decomposition
	"GeneralizedSchur"	based on generalized Schur decomposition
	"InverseFree"	a variant of "GeneralizedSchur"
	"MatrixSign"	iterative method using the matrix sign function
	"Newton"	iterative Newton method
	"Schur"	based on Schur decomposition

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

Examples

open all close all

Basic Examples (1)

Solve a discrete algebraic Riccati equation:

Wolfram Language code:

{a, b, q, r} = {(⁠|   |    |   |
| - | -- | - |
| 1 | -1 | 1 |
| 0 | 1  | 1 |
| 0 | 0  | 1 |⁠), (⁠|   |   |
| - | - |
| 1 | 0 |
| 1 | 0 |
| 0 | 1 |⁠), (⁠|    |   |     |
| -- | - | --- |
| 10 | 0 | 0   |
| 0  | 1 | 0   |
| 0  | 0 | 0.1 |⁠), (⁠|    |     |
| -- | --- |
| 10 | 0   |
| 0  | 0.1 |⁠)};

Wolfram Language code: (x = DiscreteRiccatiSolve[{a, b}, {q, r}])//MatrixForm

Verify the solution:

Wolfram Language code: Transpose[a].x.a - x - Transpose[a].x.b.Inverse[r + Transpose[b].x.b].Transpose[b].x.a + q//Chop

Scope (3)

Solve a discrete Riccati equation:

Wolfram Language code:

DiscreteRiccatiSolve[{(|      |     |
| ---- | --- |
| -1   | 0.5 |
| -0.5 | -2  |), (|     |
| --- |
| 1   |
| 0.5 |)}, {(|   |   |
| - | - |
| 1 | 0 |
| 0 | 2 |), {{1}}}]//MatrixForm

Solve a discrete Riccati equation with a cross-coupling matrix:

Wolfram Language code:

DiscreteRiccatiSolve[{(|   |        |     |
| - | ------ | --- |
| 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 |)}]//MatrixForm

Solve the Riccati equation for a sampled-data transfer-function model:

Wolfram Language code:

tfm = TransferFunctionModel[{{{1}}, (-1 + z)*(-0.605 + z)}, z, 
 SamplingPeriod -> 0.1];

Wolfram Language code:

DiscreteRiccatiSolve[Normal[StateSpaceModel[tfm]][[1 ;; 2]], {(|    |     |
| -- | --- |
| 10 | 0   |
| 0  | 0.1 |), (⁠1⁠)}]

Options (7)

Method (7)

Automatic and "Eigensystem" methods can be used for exact systems:

Wolfram Language code:

x1 = DiscreteRiccatiSolve[{(|   |   |
| - | - |
| 0 | 2 |
| 1 | 0 |), (|   |
| - |
| 0 |
| 1 |)}, {(|   |       |
| - | ----- |
| 1 | 0     |
| 0 | (4/9) |), (⁠1⁠)}, Method -> Automatic]

Wolfram Language code:

x2 = DiscreteRiccatiSolve[{(|   |   |
| - | - |
| 0 | 2 |
| 1 | 0 |), (|   |
| - |
| 0 |
| 1 |)}, {(|   |       |
| - | ----- |
| 1 | 0     |
| 0 | (4/9) |), (⁠1⁠)}, Method -> "Eigensystem"]

Wolfram Language code: x1 == x2

As well as for inexact systems:

Wolfram Language code:

x1 = DiscreteRiccatiSolve[N@{(|   |   |
| - | - |
| 0 | 2 |
| 1 | 0 |), (|   |
| - |
| 0 |
| 1 |)}, {(|   |       |
| - | ----- |
| 1 | 0     |
| 0 | (4/9) |), (⁠1⁠)}, Method -> Automatic]//Chop

Wolfram Language code:

x2 = DiscreteRiccatiSolve[N@{(|   |   |
| - | - |
| 0 | 2 |
| 1 | 0 |), (|   |
| - |
| 0 |
| 1 |)}, {(|   |       |
| - | ----- |
| 1 | 0     |
| 0 | (4/9) |), (⁠1⁠)}, Method -> "Eigensystem"]//Chop

Wolfram Language code: x1 == x2

"Schur" method can be used for inexact systems:

Wolfram Language code:

DiscreteRiccatiSolve[{(|   |   |
| - | - |
| 0 | 2 |
| 1 | 1 |), (|   |
| - |
| 0 |
| 1 |)}, {(|      |      |
| ---- | ---- |
| 2.25 | -1.5 |
| -1.5 | 1    |), (⁠2⁠)}, Method -> "Schur"]

"Newton" applies to inexact systems and may be more accurate than Automatic:

Wolfram Language code:

a = SparseArray[{{1, 1} -> 0.8, {1, 5} -> -0.8, {2, 1} -> 0.8, {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:

Wolfram Language code: x1 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "Newton"];

Wolfram Language code: Norm[a^.x1.a - x1 - a^.x1.b.Inverse[r + b^.x1.b].b^.x1.a + q]

Compare with the default method:

Wolfram Language code: x2 = DiscreteRiccatiSolve[{a, b}, {q, r}];

Wolfram Language code: Norm[a^.x2.a - x2 - a^.x2.b.Inverse[r + b^.x2.b].b^.x2.a + q]

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

Wolfram Language code:

x3 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> {"Newton", "StartingMatrix" -> x2, "Tolerance" -> 10^-15, "MaxIterations" -> 10}];

Wolfram Language code: Norm[a^.x3.a - x3 - a^.x3.b.Inverse[r + b^.x3.b].b^.x3.a + q]

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

Wolfram Language code:

a = {{0., -1.}, {0., 2.}};
b = {{1, 0}, {1, 1}};
q = {{1, 0}, {0, 0}};
r = {{4, 2}, {2, 1}};
x1 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "Newton"];

Compare with the Automatic method:

Wolfram Language code:

x2 = DiscreteRiccatiSolve[{a, b}, {q, r}];
Norm[a^.x2.a - x2 - a^.x2.b.Inverse[r + b^.x2.b].b^.x2.a + q]

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

Wolfram Language code:

{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 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "MatrixSign"];

Use xinit to initialize the "Newton" method:

Wolfram Language code: x = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> {"Newton", "StartingMatrix" -> xinit, "Tolerance" -> 10^-12}];

Compare errors:

Wolfram Language code: Map[Norm[a^.#.a - # - a^.#.b.Inverse[r + b^.#.b].b^.#.a + q]&, {xinit, x}]

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

Wolfram Language code: x3 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> {"MatrixSign", "Tolerance" -> 10^-12, "MaxIterations" -> 10}];

Wolfram Language code: Norm[a^.x3.a - x3 - a^.x3.b.Inverse[r + b^.x3.b].b^.x3.a + q]

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

Wolfram Language code:

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

The matrix a is singular:

Wolfram Language code: Det[a]

These two methods work for a singular a:

Wolfram Language code: x1 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "GeneralizedSchur"];

Wolfram Language code: x2 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "GeneralizedEigensystem"];

Verify error:

Wolfram Language code: Map[Norm[a^.#.a - # - a^.#.b.Inverse[r + b^.#.b].b^.#.a + q]&, {x1, x2}]

"Eigensystem" cannot be used when a is singular:

Wolfram Language code: x3 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "Eigensystem"];

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

Wolfram Language code:

{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.}}};

The matrix r has a high condition number:

Wolfram Language code: Norm[r] Norm[Inverse[r]]

Use the method:

Wolfram Language code: x1 = DiscreteRiccatiSolve[{a, b}, {q, r}, Method -> "InverseFree"]

Compare to the default method:

Wolfram Language code: x2 = DiscreteRiccatiSolve[{a, b}, {q, r}]

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

Wolfram Language code:

{Norm[a^.x1.a - x1 - a^.x1.b.Inverse[r + b^.x1.b].b^.x1.a + q], 
Norm[a^.x2.a - x2 - a^.x2.b.Inverse[r + b^.x2.b].b^.x2.a + q]}

Applications (2)

The linearized model of a turbo-generator has state and input matrices:

Wolfram Language code:

{a, b} = {(⁠|         |         |         |         |         |         |
| ------- | ------- | ------- | ------- | ------- | ------- |
| 0.1346  | 0.1236  | -0.0361 | 0.0037  | 0.0004  | -0.0003 |
| -0.1091 | 0.5412  | 0.3851  | -0.0631 | -0.052  | 0.0152  |
| 0.0426  | 0.1052  | 0.7915  | 0.07    | 0.0504  | -0.0172 |
| -0.0045 | 0.0205  | -0.0542 | 0.7932  | -0.5678 | 0.0025  |
| 0.0022  | -0.015  | 0.0384  | 0.5681  | 0.7526  | 0.0357  |
| 0.0002  | -0.0162 | 0.0142  | 0.0004  | -0.0299 | 0.9784  |⁠), (⁠|         |         |
| ------- | ------- |
| -0.0121 | 0.0117  |
| -0.0046 | -0.496  |
| -0.015  | 0.5517  |
| 0.0095  | -0.1763 |
| -0.0055 | -1.0216 |
| 0.0025  | 1.3944  |⁠)};

The discrete LQ regulation cost J_opt for given weighting matrices {q,r} and initial state s₀:

Wolfram Language code:

{q, r} = {(⁠|   |   |   |   |   |   |
| - | - | - | - | - | - |
| 1 | 0 | 0 | 0 | 0 | 0 |
| 0 | 1 | 0 | 0 | 0 | 0 |
| 0 | 0 | 1 | 0 | 0 | 0 |
| 0 | 0 | 0 | 1 | 0 | 0 |
| 0 | 0 | 0 | 0 | 1 | 0 |
| 0 | 0 | 0 | 0 | 0 | 1 |⁠), (⁠|   |   |
| - | - |
| 1 | 0 |
| 0 | 1 |⁠)};

Wolfram Language code: Subscript[s, 0] = {{0.1, 0, 0, 0, 0, 0}};

Wolfram Language code: x = DiscreteRiccatiSolve[{a, b}, {q, r}];

Wolfram Language code: Subscript[J, opt] = ( Subscript[s, 0].x.Transpose[Subscript[s, 0]])[[1, 1]]

Compute an optimal state-feedback gain that guarantees that all the closed-loop poles are inside a circle of radius :

Wolfram Language code:

{a, b, q, r} = {(⁠|   |          |            |
| - | -------- | ---------- |
| 1 | 0.095162 | 0.00555872 |
| 0 | 0.904837 | 0.117397   |
| 0 | 0        | 1.49182    |⁠), (⁠|             |
| ----------- |
| 0.000180325 |
| 0.00555872  |
| 0.122956    |⁠), (⁠|    |   |   |
| -- | - | - |
| 10 | 0 | 0 |
| 0  | 1 | 0 |
| 0  | 0 | 1 |⁠), (⁠1⁠)};

Wolfram Language code: {Subscript[a, 1], Subscript[b, 1]} = ({a, b}/α) /. α -> 0.7;

Wolfram Language code: x = DiscreteRiccatiSolve[{Subscript[a, 1], Subscript[b, 1]}, {q, r}];

Wolfram Language code:

Eigenvalues[a - b.Inverse[Transpose[Subscript[b, 1]].x.Subscript[b, 1] + r].Transpose[Subscript[b, 1]].x.Subscript[a, 1]]

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

Wolfram Language code: x = DiscreteRiccatiSolve[{a, b}, {q, r}];

Wolfram Language code: Eigenvalues[a - b.Inverse[Transpose[b].x.b + r].Transpose[b].x.a]

Properties & Relations (11)

DiscreteRiccatiSolve[{a,b},{q,r,p}] is equivalent to DiscreteRiccatiSolve[{a-b..p,b},{q-p.r^-1.p,r}]:

Wolfram Language code:

{a, b, q, r, p} = {(⁠|                 |                   |                |
| --------------- | ----------------- | -------------- |
| -0.25 - 0.011 I | 0.0947  + 0.089 I | 0.35  + 0.16 I |
| 0.35  + 0.22 I  | -0.08 - 0.22 I    | -0.4 - 0.53 I  |
| -0.16 - 0.07 I  | 0.05  + 0.1 I     | 0.22  + 0.05 I |⁠), (⁠|           |             |
| --------- | ----------- |
| 1 + 0.5 I | -1 - 0.5 I  |
| -0.5 - 1I | -0.5 + 1.5I |
| 0.1 I     | -0.5I       |⁠), (⁠|   |   |   |
| - | - | - |
| 1 | 0 | 0 |
| 0 | 1 | 0 |
| 0 | 0 | 1 |⁠), (⁠|   |    |
| - | -- |
| 3 | -I |
| I | -3 |⁠), (⁠|      |     |
| ---- | --- |
| 0.1  | 0   |
| -0.2 | 0   |
| 0.15 | 0.2 |⁠)};

Wolfram Language code: x = DiscreteRiccatiSolve[{a, b}, {q, r, p}];

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

Wolfram Language code:

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

Wolfram Language code: Chop[x - xp]

If {a,b} is stabilizable and {a,g} is detectable, and q=Transpose[g].g, then the solution to the discrete Riccati equation is positive semidefinite:

Wolfram Language code:

{a, b, g, r} = {(⁠|      |      |     |
| ---- | ---- | --- |
| 0.45 | 0    | 0   |
| 1    | -0.9 | 1   |
| 0    | 0    | 0.7 |⁠), (⁠|   |
| - |
| 1 |
| 1 |
| 0 |⁠), (⁠1	0	0⁠), (⁠1⁠)};

Wolfram Language code: JordanModelDecomposition[StateSpaceModel[{a, b, g}]]//Last

Wolfram Language code: PositiveSemidefiniteMatrixQ@DiscreteRiccatiSolve[{a, b}, {Transpose[g].g, r}]

If {a,b} is controllable and {a,g} is observable, and q=Transpose[g].g, then the solution to the discrete Riccati equation is positive definite:

Wolfram Language code:

{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⁠)};

Wolfram Language code: ControllableModelQ[StateSpaceModel[{a, b}]]

Wolfram Language code: ControllableModelQ[StateSpaceModel[{Transpose[a], Transpose[g]}]]

Wolfram Language code: PositiveDefiniteMatrixQ[DiscreteRiccatiSolve[{a, b}, {Transpose[g].g, r}]]

The matrix associated with the discrete algebraic Riccati equation is symplectic:

Wolfram Language code:

With[{a = (⁠|   |     |         |
| - | --- | ------- |
| 1 | 0.5 | 0.4     |
| 0 | 1   | 0.09932 |
| 0 | 0   | 0.6738  |⁠), 
	b = (⁠|         |
| ------- |
| 0.0085  |
| 0.04    |
| 0.09932 |⁠), 
	q = (⁠|   |   |   |
| - | - | - |
| 5 | 0 | 0 |
| 0 | 3 | 0 |
| 0 | 0 | 1 |⁠), r = (⁠1⁠)}, 
	M = (|                                    |                               |
| ---------------------------------- | ----------------------------- |
| a + b.Inverse[r].b.Inverse[a^].q | -b.Inverse[r].b.Inverse[a^] |
| -Inverse[a^].q                    | Inverse[a^]                  |)//ArrayFlatten];

Wolfram Language code:

Ω = ArrayFlatten[(|                          |                          |
| ------------------------ | ------------------------ |
| ConstantArray[0, {3, 3}] | IdentityMatrix[3]        |
| -IdentityMatrix[3]       | ConstantArray[0, {3, 3}] |)];

Wolfram Language code: Chop[ConjugateTranspose[M].Ω.M] == Ω

The eigenvalues of the symplectic matrix are pairs of the form :

Wolfram Language code:

With[{a = (⁠|                    |                    |                    |
| ------------------ | ------------------ | ------------------ |
| 0.0661 - 0.2022 I  | -0.8839 - 0.4677 I | -0.9869 - 0.2443 I |
| -0.7635 - 1.0728 I | 0.312  - 1.66 I    | -0.8182 - 1.546 I  |
| 1.311 + 1.3799 I   | 1.0445  + 2.1423 I | 2.2958 + 1.8621 I  |⁠), 
	b = (⁠|                    |
| ------------------ |
| 0.0015  + 0.0062 I |
| 0.0098  + 0.01223I |
| -0.004 - 0.01371I  |⁠), 
	q = (⁠|   |   |   |
| - | - | - |
| 5 | 0 | 0 |
| 0 | 3 | 0 |
| 0 | 0 | 1 |⁠), r = (⁠1⁠)}, 
	M = (|                                   |                              |
| --------------------------------- | ---------------------------- |
| a + b.Inverse[r].b.Inverse[a].q | -b.Inverse[r].b.Inverse[a] |
| -Inverse[a].q                    | Inverse[a]                  |)//ArrayFlatten];

Wolfram Language code:

ListPlot[Eigenvalues[M] /. Complex[x_, y_] :> {x, y}, PlotStyle -> PointSize[Medium], Epilog -> {Red, Circle[{0, 0}]}, AxesOrigin -> {0, 0}]

and are similar matrices:

Wolfram Language code: Complement[Eigenvalues[M], Eigenvalues[Inverse[M]], SameTest -> (Chop[#1 - #2] == 0&)]

The symplectic matrix must satisfy the stability and complementarity properties to obtain a stabilizing solution to the Riccati equation:

Wolfram Language code:

a = b = q = r = (⁠|   |   |
| - | - |
| 1 | 0 |
| 0 | 1 |⁠);

Wolfram Language code:

{vals, vecs} = Eigensystem[(|                                   |                               |
| --------------------------------- | ----------------------------- |
| a + b.Inverse[r].b.Inverse[a].q | -b.Inverse[ r].b.Inverse[a] |
| -Inverse[a].q                    | Inverse[a]                   |)//ArrayFlatten]

Stability property:

Wolfram Language code: Cases[vals, _ ? (Abs[#] == 1&)]

Complementarity property:

Wolfram Language code: stableBasis = Extract[vecs, Position[vals, _ ? (Abs[#] < 1&), {1}]];

Wolfram Language code: {{Subscript[x, 1], Subscript[x, 2]}} = Partition[stableBasis, {2, 2}]

Wolfram Language code: MatrixRank[Subscript[x, 1]] == 2

The stabilizing solution:

Wolfram Language code: Subscript[x, 2].Inverse[Subscript[x, 1]]//Simplify

Wolfram Language code: DiscreteRiccatiSolve[{a, b}, {q, r}] == %

The eigenvalues of the feedback system with are the stable eigenvalues of the symplectic matrix:

Wolfram Language code:

{a, b, q, r} = {(⁠|   |   |
| - | - |
| 1 | 3 |
| 1 | 2 |⁠), (⁠|    |
| -- |
| 1  |
| -2 |⁠), (⁠|     |   |
| --- | - |
| 0.2 | 0 |
| 0   | 1 |⁠), (⁠0.5)};

Wolfram Language code: With[{x = DiscreteRiccatiSolve[{a, b}, {q, r}]}, Eigenvalues[a - b.Inverse[Transpose[b].x.b + r].Transpose[b].x.a]]

Wolfram Language code:

Eigenvalues[ArrayFlatten[(|                                   |                               |
| --------------------------------- | ----------------------------- |
| a + b.Inverse[r].b.Inverse[a].q | -b.Inverse[ r].b.Inverse[a] |
| -Inverse[a].q                    | Inverse[a]                   |)]];
Select[%, Abs[#] < 1&]

Compute optimal state-feedback gains using DiscreteRiccatiSolve:

Wolfram Language code:

{a, b, q, r} = {(⁠|     |     |
| --- | --- |
| 0.2 | 1   |
| -1  | 0.5 |⁠), (⁠|   |
| - |
| 1 |
| 0 |⁠), (⁠|   |   |
| - | - |
| 1 | 0 |
| 0 | 1 |⁠), (⁠1⁠)};

Wolfram Language code: With[{x = DiscreteRiccatiSolve[{a, b}, {q, r}]}, Inverse[Transpose[b].x.b + r].Transpose[b].x.a]//N

Use LQRegulatorGains to obtain the same result directly:

Wolfram Language code: LQRegulatorGains[StateSpaceModel[{a, b}, SamplingPeriod -> T], {q, r}]//N

Compute optimal output feedback gains using DiscreteRiccatiSolve:

Wolfram Language code:

{a, b, c, q, r} = {(⁠|    |    |
| -- | -- |
| -1 | 1  |
| 0  | -1 |⁠), (⁠|   |   |
| - | - |
| 0 | 1 |
| 1 | 0 |⁠), (⁠1	0⁠), (⁠10⁠), (⁠|   |   |
| - | - |
| 1 | 0 |
| 0 | 1 |⁠)};

Wolfram Language code:

With[{x = DiscreteRiccatiSolve[{a, b}, {Transpose[c].q.c, r}]}, Inverse[Transpose[b].x.b + r].Transpose[b].x.a]//N//MatrixForm

LQOutputRegulatorGains gives the same result:

Wolfram Language code: LQOutputRegulatorGains[StateSpaceModel[{a, b, c}, SamplingPeriod -> T], {q, r}]//N//MatrixForm

Compute optimal estimator gains using DiscreteRiccatiSolve:

Wolfram Language code:

{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)};

Wolfram Language code:

{w, v} = {(⁠|       |       |
| ----- | ----- |
| 10^-2 | 0     |
| 0     | 10^-2 |⁠), (⁠10^-4⁠)};

Wolfram Language code:

With[{x = DiscreteRiccatiSolve[{Transpose[a], Transpose[c]}, {b.w.Transpose[b], v}]}, x.Transpose[c].Inverse[c.x.Transpose[c] + v]]

Use LQEstimatorGains:

Wolfram Language code: LQEstimatorGains[StateSpaceModel[{a, b, c}, SamplingPeriod -> τ], {w, v}]

The optimal trajectory of the discrete approximation of a system results in a higher cost:

Wolfram Language code:

ssm = StateSpaceModel[{{{0, 0, 1, 0}, {0, 0, 0, 1}, {10.7275, 0, 0, 0}, {-0.9275, 0, 0, 0}}, 
  {{0}, {0}, {-0.178571}, {0.178571}}, {{1, 0, 0, 0}, {0, 1, 0, 0}}, {{0}, {0}}}, 
 SamplingPeriod -> None, SystemsModelLabels -> None];

Wolfram Language code: ssmd = ToDiscreteTimeModel[ssm, 0.1];

Wolfram Language code:

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

The optimal LQ regulation costs for the systems starting at s₀:

Wolfram Language code: Subscript[s, 0] = {{0}, {0.1}, {0}, {0}};

Wolfram Language code: J = (Transpose[Subscript[s, 0]].RiccatiSolve[Normal[ssm][[1 ;; 2]], {q, r}].Subscript[s, 0])[[1, 1]];

Wolfram Language code:

Subscript[J, d] = (Transpose[Subscript[s, 0]].DiscreteRiccatiSolve[Normal[ssmd][[1 ;; 2]], {q, r}].Subscript[s, 0])[[1, 1]];

Wolfram Language code: Greater[Subscript[J, d], J]

Possible Issues (1)

If {a,b} is not stabilizable and {a,g} is not detectable, then the Riccati equation with q=g.g has no stabilizing solution:

Wolfram Language code:

{a, b, q, g, r} = {(⁠|     |     |      |
| --- | --- | ---- |
| 0.3 | 0   | 0    |
| 0.3 | 0.4 | 0.06 |
| 0   | 0   | 1.2  |⁠), (⁠|       |
| ----- |
| 0.006 |
| 0.04  |
| 0     |⁠), (⁠|   |   |   |
| - | - | - |
| 1 | 0 | 0 |
| 0 | 0 | 0 |
| 0 | 0 | 0 |⁠), (⁠1	0	0⁠), (⁠1⁠)};

Wolfram Language code: {ControllableModelQ[StateSpaceModel[{a, b}]], ControllableModelQ[StateSpaceModel[{Transpose[a], Transpose[g]}]]}

Wolfram Language code: DiscreteRiccatiSolve[{a, b}, {q, r}]

Top

More Learning

Tech Support

Wolfram Solutions

Wolfram Solutions For Education

Get Started

Grow Your Skills

Work with Us

Educational Programs for Adults

Educational Programs for Youth

Read

DiscreteRiccatiSolve

Details and Options

Examples

Basic Examples (1)

Scope (3)

Options (7)

Method (7)

Applications (2)

Properties & Relations (11)

Possible Issues (1)

Text

CMS

APA

BibTeX

BibLaTeX

DiscreteRiccatiSolve

Details and Options

Examples

Basic Examples (1)

Scope (3)

Options (7)

Method (7)

Applications (2)

Properties & Relations (11)

Possible Issues (1)

See Also

Related Guides

History

Text

CMS

APA

BibTeX

BibLaTeX