Skip to content

Commit 05a419d

Browse files
tamaskistamaskis
committed
Updated matrix ODE examples and documentation
1 parent 117dfc3 commit 05a419d

24 files changed

+74
-58
lines changed

‎docs/Matrix_ODE_Example_doc.html‎

Lines changed: 29 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
<!--
77
This HTML was auto-generated from MATLAB code.
88
To make changes, update the MATLAB code and republish this document.
9-
--><title>Matrix ODE Example</title><meta name="generator" content="MATLAB 9.11"><link rel="schema.DC" href="http://purl.org/dc/elements/1.1/"><meta name="DC.date" content="2021-12-22"><meta name="DC.source" content="Matrix_ODE_Example_doc.m"><style type="text/css">
9+
--><title>Matrix ODE Example</title><meta name="generator" content="MATLAB 9.11"><link rel="schema.DC" href="http://purl.org/dc/elements/1.1/"><meta name="DC.date" content="2021-12-23"><meta name="DC.source" content="Matrix_ODE_Example_doc.m"><style type="text/css">
1010
html,body,div,span,applet,object,iframe,h1,h2,h3,h4,h5,h6,p,blockquote,pre,a,abbr,acronym,address,big,cite,code,del,dfn,em,font,img,ins,kbd,q,s,samp,small,strike,strong,tt,var,b,u,i,center,dl,dt,dd,ol,ul,li,fieldset,form,label,legend,table,caption,tbody,tfoot,thead,tr,th,td{margin:0;padding:0;border:0;outline:0;font-size:100%;vertical-align:baseline;background:transparent}body{line-height:1}ol,ul{list-style:none}blockquote,q{quotes:none}blockquote:before,blockquote:after,q:before,q:after{content:'';content:none}:focus{outine:0}ins{text-decoration:none}del{text-decoration:line-through}table{border-collapse:collapse;border-spacing:0}
1111

1212
html { min-height:100%; margin-bottom:1px; }
@@ -67,42 +67,41 @@
6767

6868

6969

70-
</style></head><body><div class="content"><h1>Matrix ODE Example</h1><!--introduction--><p><a href="index.html">Back to ODE Solver Toolbox Contents</a>.</p><!--/introduction--><h2>Contents</h2><div><ul><li><a href="#1">Problem Statement</a></li><li><a href="#2">Matrices</a></li><li><a href="#3">Final Condition</a></li><li><a href="#4"><tt>odefun_mat2vec</tt></a></li><li><a href="#7"><tt>odeIC_mat2vec</tt></a></li><li><a href="#9"><tt>odesol_vec2mat</tt></a></li><li><a href="#12">Solution for <img src="Matrix_ODE_Example_doc_eq01737430298148732982.png" alt="$\mathbf{P}_{0}$" style="width:12px;height:10px;"></a></li><li><a href="#13">References</a></li></ul></div><h2 id="1">Problem Statement</h2><p>Consider the Riccati differential equation for the finite-horizon linear quadratic regulator problem:</p><p><img src="Matrix_ODE_Example_doc_eq18025082327632211741.png" alt="$$\dot{\mathbf{P}}=-\left[\mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}-(\mathbf{P}\mathbf{B}+\mathbf{N})\mathbf{R}^{-1}(\mathbf{B}^{T}\mathbf{P}+\mathbf{N}^{T})+\mathbf{Q}\right]$$" style="width:263px;height:15px;"></p><p>Find <img src="Matrix_ODE_Example_doc_eq11156567409250041374.png" alt="$\mathbf{P}(0)=\mathbf{P}_{0}$" style="width:50px;height:11px;"> given that</p><p><img src="Matrix_ODE_Example_doc_eq06697316367777428526.png" alt="$$\mathbf{P}(T)=\mathbf{P}_{T}=\pmatrix{1&amp;1\cr1&amp;1}$$" style="width:109px;height:27px;"></p><p>where <img src="Matrix_ODE_Example_doc_eq04333988821981505429.png" alt="$T=5$" style="width:28px;height:8px;">. The state (<img src="Matrix_ODE_Example_doc_eq15190527482521674856.png" alt="$\mathbf{A}$" style="width:9px;height:8px;">) and input (<img src="Matrix_ODE_Example_doc_eq16427099040942302163.png" alt="$\mathbf{B}$" style="width:9px;height:8px;">) matrices are</p><p><img src="Matrix_ODE_Example_doc_eq08594218919394792437.png" alt="$$\mathbf{A}=\pmatrix{1&amp;1\cr2&amp;1},\quad\quad\mathbf{B}=\pmatrix{1\cr1}$$" style="width:142px;height:27px;"></p><p>The cross-coupling (<img src="Matrix_ODE_Example_doc_eq15866111403695682364.png" alt="$\mathbf{N}$" style="width:10px;height:8px;">), state (<img src="Matrix_ODE_Example_doc_eq04980437671470961831.png" alt="$\mathbf{Q}$" style="width:9px;height:10px;">), and input (<img src="Matrix_ODE_Example_doc_eq07378229934469313194.png" alt="$\mathbf{R}$" style="width:10px;height:8px;">) weighting matrices are</p><p><img src="Matrix_ODE_Example_doc_eq04665828541583697595.png" alt="$$\mathbf{N}=\pmatrix{0&amp;0\cr0&amp;0},\quad\quad\mathbf{Q}=\pmatrix{2&amp;1\cr1&amp;1},\quad\quad\mathbf{R}=\pmatrix{1&amp;0\cr0&amp;1}$$" style="width:254px;height:27px;"></p><h2 id="2">Matrices</h2><pre class="codeinput"><span class="comment">% state matrix</span>
70+
</style></head><body><div class="content"><h1>Matrix ODE Example</h1><!--introduction--><p><a href="index.html">Back to ODE Solver Toolbox Contents</a>.</p><!--/introduction--><h2>Contents</h2><div><ul><li><a href="#1">Problem Statement</a></li><li><a href="#2">Matrices</a></li><li><a href="#3">Terminal Condition</a></li><li><a href="#4"><tt>odefun_mat2vec</tt></a></li><li><a href="#7"><tt>odeIC_mat2vec</tt></a></li><li><a href="#9"><tt>odesol_vec2mat</tt></a></li><li><a href="#12">Solution for <img src="Matrix_ODE_Example_doc_eq01737430298148732982.png" alt="$\mathbf{P}_{0}$" style="width:12px;height:10px;"></a></li><li><a href="#13">References</a></li></ul></div><h2 id="1">Problem Statement</h2><p>Consider the Riccati differential equation for the finite-horizon linear quadratic regulator problem:</p><p><img src="Matrix_ODE_Example_doc_eq15499043028313771354.png" alt="$$\dot{\mathbf{P}}=-\left[\mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}-(\mathbf{P}\mathbf{B}+\mathbf{S})\mathbf{R}^{-1}(\mathbf{B}^{T}\mathbf{P}+\mathbf{S}^{T})+\mathbf{Q}\right]$$" style="width:258px;height:15px;"></p><p>Find <img src="Matrix_ODE_Example_doc_eq11156567409250041374.png" alt="$\mathbf{P}(0)=\mathbf{P}_{0}$" style="width:50px;height:11px;"> given that</p><p><img src="Matrix_ODE_Example_doc_eq06697316367777428526.png" alt="$$\mathbf{P}(T)=\mathbf{P}_{T}=\pmatrix{1&amp;1\cr1&amp;1}$$" style="width:109px;height:27px;"></p><p>where <img src="Matrix_ODE_Example_doc_eq04333988821981505429.png" alt="$T=5$" style="width:28px;height:8px;">. The state (<img src="Matrix_ODE_Example_doc_eq15190527482521674856.png" alt="$\mathbf{A}$" style="width:9px;height:8px;">) and input (<img src="Matrix_ODE_Example_doc_eq16427099040942302163.png" alt="$\mathbf{B}$" style="width:9px;height:8px;">) matrices are</p><p><img src="Matrix_ODE_Example_doc_eq08594218919394792437.png" alt="$$\mathbf{A}=\pmatrix{1&amp;1\cr2&amp;1},\quad\quad\mathbf{B}=\pmatrix{1\cr1}$$" style="width:142px;height:27px;"></p><p>The state (<img src="Matrix_ODE_Example_doc_eq04980437671470961831.png" alt="$\mathbf{Q}$" style="width:9px;height:10px;">), input (<img src="Matrix_ODE_Example_doc_eq07378229934469313194.png" alt="$\mathbf{R}$" style="width:10px;height:8px;">), and cross-coupling (<img src="Matrix_ODE_Example_doc_eq17705281529142768610.png" alt="$\mathbf{S}$" style="width:6px;height:8px;">) weighting matrices are</p><p><img src="Matrix_ODE_Example_doc_eq06656902451109221026.png" alt="$$\mathbf{Q}=\pmatrix{2&amp;1\cr1&amp;1},\quad\quad\mathbf{R}=\pmatrix{1},\quad\quad\mathbf{S}=\pmatrix{0\cr0}$$" style="width:210px;height:27px;"></p><h2 id="2">Matrices</h2><pre class="codeinput"><span class="comment">% state matrix</span>
7171
A = [1 1;
7272
2 1];
7373

7474
<span class="comment">% input matrix</span>
7575
B = [1;
7676
1];
7777

78-
<span class="comment">% cross-coupling weighting matrix</span>
79-
N = [0 0;
80-
0 0];
81-
8278
<span class="comment">% state weighting matrix</span>
8379
Q = [2 1;
8480
1 1];
8581

8682
<span class="comment">% input weighting matrix</span>
87-
R = [1 0;
88-
0 1];
89-
</pre><h2 id="3">Final Condition</h2><pre class="codeinput"><span class="comment">% final condition</span>
83+
R = 1;
84+
85+
<span class="comment">% cross-coupling weighting matrix</span>
86+
S = [0;
87+
0];
88+
</pre><h2 id="3">Terminal Condition</h2><pre class="codeinput"><span class="comment">% terminal condition</span>
9089
PT = [1 1;
9190
1 1];
9291

9392
<span class="comment">% final time</span>
9493
T = 5;
95-
</pre><h2 id="4"><tt>odefun_mat2vec</tt></h2><p>Defining the Riccati differential equation (a matrix-valued ODE),</p><pre class="codeinput">F = @(t,P) -(A.'*P+P*A-(P*B+N)/R*(B.'*P+N.')+Q);
94+
</pre><h2 id="4"><tt>odefun_mat2vec</tt></h2><p>Defining the Riccati differential equation (a matrix-valued ODE),</p><pre class="codeinput">F = @(t,P) -(A.'*P+P*A-(P*B+S)/R*(B.'*P+S.')+Q);
9695
</pre><p>Converting this matrix-valued ODE to a vector-valued ODE,</p><pre class="codeinput">f = odefun_mat2vec(F);
97-
</pre><p>Note that since <img src="Matrix_ODE_Example_doc_eq01284226214242475479.png" alt="$\mathbf{P}$" style="width:8px;height:8px;"> is a square matrix, we did not have to specify its number of rows for the <tt>odefun_mat2vec</tt> function.</p><h2 id="7"><tt>odeIC_mat2vec</tt></h2><p>In this case, we have a final condition at time <img src="Matrix_ODE_Example_doc_eq16373138058832967206.png" alt="$t=T=5$" style="width:47px;height:8px;">, and want to find the initial condition at time <img src="Matrix_ODE_Example_doc_eq02656990568186096865.png" alt="$t=t_{0}=0$" style="width:47px;height:10px;">. To do this, we will integrate backwards using an ODE solver, so the final condition <img src="Matrix_ODE_Example_doc_eq07282847098330546007.png" alt="$\mathbf{P}_{T}$" style="width:14px;height:10px;"> will actually form the initial condition. Therefore, finding the final condition for the corresponding vector-valued ODE,</p><pre class="codeinput">yT = odeIC_mat2vec(PT);
96+
</pre><p>Note that since <img src="Matrix_ODE_Example_doc_eq01284226214242475479.png" alt="$\mathbf{P}$" style="width:8px;height:8px;"> is a square matrix, we did not have to specify its number of rows for the <tt>odefun_mat2vec</tt> function.</p><h2 id="7"><tt>odeIC_mat2vec</tt></h2><p>In this case, we have a terminal condition at time <img src="Matrix_ODE_Example_doc_eq16373138058832967206.png" alt="$t=T=5$" style="width:47px;height:8px;">, and want to find the initial condition at time <img src="Matrix_ODE_Example_doc_eq02656990568186096865.png" alt="$t=t_{0}=0$" style="width:47px;height:10px;">. To do this, we will integrate backwards using an ODE solver, so the terminal condition <img src="Matrix_ODE_Example_doc_eq07282847098330546007.png" alt="$\mathbf{P}_{T}$" style="width:14px;height:10px;"> will actually form the initial condition. Therefore, finding the terminal condition for the corresponding vector-valued ODE,</p><pre class="codeinput">yT = odeIC_mat2vec(PT);
9897
</pre><p>Note that since <img src="Matrix_ODE_Example_doc_eq01284226214242475479.png" alt="$\mathbf{P}$" style="width:8px;height:8px;"> is a square matrix, we did not have to specify its number of rows for the <tt>odeIC_mat2vec</tt> function.</p><h2 id="9"><tt>odesol_vec2mat</tt></h2><p>First, let's solve the vector-valued ODE for <img src="Matrix_ODE_Example_doc_eq13802958261507554970.png" alt="$\mathbf{y}(t)$" style="width:19px;height:11px;">. Solving the ODE from <img src="Matrix_ODE_Example_doc_eq12231604572696234117.png" alt="$t=T$" style="width:27px;height:8px;"> to <img src="Matrix_ODE_Example_doc_eq00402237798869806649.png" alt="$t=0$" style="width:24px;height:8px;">,</p><pre class="codeinput">[t,y] = ode45(f,[T,0],yT);
9998
</pre><p>Transforming the solution matrix (<tt>y</tt>) for the vector-valued ODE into the solution <i>array</i> (<tt>P</tt>) for the matrix-valued ODE,</p><pre class="codeinput">P = odesol_vec2mat(y);
10099
</pre><p>Note that since <img src="Matrix_ODE_Example_doc_eq01284226214242475479.png" alt="$\mathbf{P}$" style="width:8px;height:8px;"> is a square matrix, we did not have to specify its number of rows for the <tt>odesol_vec2mat</tt> function.</p><h2 id="12">Solution for <img src="Matrix_ODE_Example_doc_eq01737430298148732982.png" alt="$\mathbf{P}_{0}$" style="width:12px;height:10px;"></h2><p>Our original goal was to find <img src="Matrix_ODE_Example_doc_eq01737430298148732982.png" alt="$\mathbf{P}_{0}$" style="width:12px;height:10px;">. Extracting this matrix from the solution array (noting that the solution corresponding to <img src="Matrix_ODE_Example_doc_eq00402237798869806649.png" alt="$t=0$" style="width:24px;height:8px;"> will be stored at the end of the solution array since we integrated backwards in time), we find</p><pre class="codeinput">P0 = P(:,:,end)
101100
</pre><pre class="codeoutput">
102101
P0 =
103102

104-
1.292500070902461 0.548255578137689
105-
0.548255578137689 0.816900870990683
103+
2.1470 1.1582
104+
1.1582 1.2546
106105

107106
</pre><h2 id="13">References</h2><p>This example is adapted from the following sources:</p><div><ul><li><a href="https://www.mathworks.com/matlabcentral/answers/94722-how-can-i-solve-the-matrix-riccati-differential-equation-within-matlab">https://www.mathworks.com/matlabcentral/answers/94722-how-can-i-solve-the-matrix-riccati-differential-equation-within-matlab</a></li><li><a href="https://en.wikipedia.org/wiki/Linear-quadratic_regulator">https://en.wikipedia.org/wiki/Linear-quadratic_regulator</a></li></ul></div><p class="footer"><br><a href="https://www.mathworks.com/products/matlab/">Published with MATLAB&reg; R2021b</a><br></p></div><!--
108107
##### SOURCE BEGIN #####
@@ -113,7 +112,7 @@
113112
% Consider the Riccati differential equation for the finite-horizon linear
114113
% quadratic regulator problem:
115114
%
116-
% $$\dot{\mathbf{P}}=-\left[\mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}-(\mathbf{P}\mathbf{B}+\mathbf{N})\mathbf{R}^{-1}(\mathbf{B}^{T}\mathbf{P}+\mathbf{N}^{T})+\mathbf{Q}\right]$$
115+
% $$\dot{\mathbf{P}}=-\left[\mathbf{A}^{T}\mathbf{P}+\mathbf{P}\mathbf{A}-(\mathbf{P}\mathbf{B}+\mathbf{S})\mathbf{R}^{-1}(\mathbf{B}^{T}\mathbf{P}+\mathbf{S}^{T})+\mathbf{Q}\right]$$
117116
%
118117
% Find $\mathbf{P}(0)=\mathbf{P}_{0}$ given that
119118
%
@@ -124,10 +123,10 @@
124123
%
125124
% $$\mathbf{A}=\pmatrix{1&1\cr2&1},\quad\quad\mathbf{B}=\pmatrix{1\cr1}$$
126125
%
127-
% The cross-coupling ($\mathbf{N}$), state ($\mathbf{Q}$), and input
128-
% ($\mathbf{R}$) weighting matrices are
126+
% The state ($\mathbf{Q}$), input ($\mathbf{R}$), and cross-coupling
127+
% ($\mathbf{S}$) weighting matrices are
129128
%
130-
% $$\mathbf{N}=\pmatrix{0&0\cr0&0},\quad\quad\mathbf{Q}=\pmatrix{2&1\cr1&1},\quad\quad\mathbf{R}=\pmatrix{1&0\cr0&1}$$
129+
% $$\mathbf{Q}=\pmatrix{2&1\cr1&1},\quad\quad\mathbf{R}=\pmatrix{1},\quad\quad\mathbf{S}=\pmatrix{0\cr0}$$
131130
132131
%% Matrices
133132
@@ -139,40 +138,39 @@
139138
B = [1;
140139
1];
141140
142-
% cross-coupling weighting matrix
143-
N = [0 0;
144-
0 0];
145-
146141
% state weighting matrix
147142
Q = [2 1;
148143
1 1];
149144
150145
% input weighting matrix
151-
R = [1 0;
152-
0 1];
153-
%% Final Condition
146+
R = 1;
147+
148+
% cross-coupling weighting matrix
149+
S = [0;
150+
0];
151+
%% Terminal Condition
154152
155-
% final condition
153+
% terminal condition
156154
PT = [1 1;
157155
1 1];
158156
159157
% final time
160158
T = 5;
161159
%% |odefun_mat2vec|
162160
% Defining the Riccati differential equation (a matrix-valued ODE),
163-
F = @(t,P) -(A.'*P+P*A-(P*B+N)/R*(B.'*P+N.')+Q);
161+
F = @(t,P) -(A.'*P+P*A-(P*B+S)/R*(B.'*P+S.')+Q);
164162
%%
165163
% Converting this matrix-valued ODE to a vector-valued ODE,
166164
f = odefun_mat2vec(F);
167165
%%
168166
% Note that since $\mathbf{P}$ is a square matrix, we did not have to
169167
% specify its number of rows for the |odefun_mat2vec| function.
170168
%% |odeIC_mat2vec|
171-
% In this case, we have a final condition at time $t=T=5$, and want to find
172-
% the initial condition at time $t=t_{0}=0$. To do this, we will integrate
173-
% backwards using an ODE solver, so the final condition $\mathbf{P}_{T}$
174-
% will actually form the initial condition. Therefore, finding the final
175-
% condition for the corresponding vector-valued ODE,
169+
% In this case, we have a terminal condition at time $t=T=5$, and want to
170+
% find the initial condition at time $t=t_{0}=0$. To do this, we will
171+
% integrate backwards using an ODE solver, so the terminal condition
172+
% $\mathbf{P}_{T}$ will actually form the initial condition. Therefore,
173+
% finding the terminal condition for the corresponding vector-valued ODE,
176174
yT = odeIC_mat2vec(PT);
177175
%%
178176
% Note that since $\mathbf{P}$ is a square matrix, we did not have to

‎docs/Matrix_ODE_Example_doc_eq00402237798869806649.png‎

0 Bytes
Loading

‎docs/Matrix_ODE_Example_doc_eq01284226214242475479.png‎

0 Bytes
Loading

‎docs/Matrix_ODE_Example_doc_eq01737430298148732982.png‎

0 Bytes
Loading

‎docs/Matrix_ODE_Example_doc_eq02656990568186096865.png‎

0 Bytes
Loading

‎docs/Matrix_ODE_Example_doc_eq04333988821981505429.png‎

0 Bytes
Loading

‎docs/Matrix_ODE_Example_doc_eq04980437671470961831.png‎

0 Bytes
Loading

‎docs/Matrix_ODE_Example_doc_eq05545325794731378589.png‎

9.68 KB
Loading

‎docs/Matrix_ODE_Example_doc_eq06656902451109221026.png‎

7.71 KB
Loading

‎docs/Matrix_ODE_Example_doc_eq06697316367777428526.png‎

0 Bytes
Loading

0 commit comments

Comments
 (0)