Before applying dynamic programming to solve the problem, the solution was first found by brute force in order to verify that the D.P. method when completed was correct. Using brute force, the optimal arrangement is found to be:
The brute force method also generated a list of all the arrangements (44 of them) and the cost of each. For reference, here is the complete table with the small Matlab code used to generated it in the appendix. After this, the graphical D.P. method called branch and bound was used to verify this result.
| district 1 | district 2 | district 3 | cost in millions |
| 0 | 0 | 0 | 4.0 |
| 0 | 0 | 1 | 3.5 |
| 0 | 0 | 2 | 3.2 |
| 0 | 0 | 3 | 2.8 |
| 0 | 1 | 0 | 3.8 |
| 0 | 1 | 1 | 3.3 |
| 0 | 1 | 2 | 3.0 |
| 0 | 1 | 3 | 2.6 |
| 0 | 2 | 0 | 3.7 |
| 0 | 2 | 1 | 3.2 |
| 0 | 2 | 2 | 2.9 |
| 0 | 2 | 3 | 2.5 |
| 0 | 3 | 0 | 3.6 |
| 0 | 3 | 1 | 3.1 |
| 0 | 3 | 2 | 2.8 |
| 1 | 0 | 0 | 2.9 |
| 1 | 0 | 1 | 2.4 |
| 1 | 0 | 2 | 2.1 |
| 1 | 0 | 3 | 1.7 |
| 1 | 1 | 0 | 2.7 |
| 1 | 1 | 1 | 2.2 |
| 1 | 1 | 2 | 1.9 |
| 1 | 1 | 3 | 1.5 |
| 1 | 2 | 0 | 2.6 |
| 1 | 2 | 1 | 2.1 |
| 1 | 2 | 2 | 1.8 |
| 1 | 3 | 0 | 2.5 |
| 1 | 3 | 1 | 2.0 |
| 2 | 0 | 0 | 2.3 |
| 2 | 0 | 1 | 1.8 |
| 2 | 0 | 2 | 1.5 |
| 2 | 0 | 3 | \(1.1^{\ast }\) |
| 2 | 1 | 0 | 2.1 |
| 2 | 1 | 1 | 1.6 |
| 2 | 1 | 2 | 1.3 |
| 2 | 2 | 0 | 2.0 |
| 2 | 2 | 1 | 1.5 |
| 2 | 3 | 0 | 1.9 |
| 3 | 0 | 0 | 2.2 |
| 3 | 0 | 1 | 1.7 |
| 3 | 0 | 2 | 1.4 |
| 3 | 1 | 0 | 2.0 |
| 3 | 1 | 1 | 1.5 |
| 3 | 2 | 0 | 1.9 |
Let the state \(x\) be the number of stations available to be assigned at each stage. For example, if we are at stage \(2\) and \(x=5\), this means all \(5\) stations are available to be assigned to district two. Stage one was the decision to decide on district one, stage two for the decision to assign for district two and the final stage, stage three is for district three. This is arbitrary, any order will give the same answer. One can decide on district three first, and then district one for example. This makes no difference to the final result.
The following diagram shows the result found which agrees with the brute force method above. The branch cost is the number above the arrow itself. The number in the small rectangle at the node, is the minimal cost of the branch leaving that node. For example, in stage three, when \(x=5\), there are 4 branches that leave that node where (0,1,2,3) stations can be assigned. The lowest cost of these branches is the one with cost \(0.3\) and that is what goes in the small square next to the node. This process continues moving backward. This method is the graphical equivalent to the Bellman dynamic equations and can be used when the number of states is finite and the number of decisions at each state is finite also.
The constraint now is
\begin{equation} 3u\left ( 0\right ) +u\left ( 1\right ) +2u\left ( 2\right ) \leq 9 \tag{1} \end{equation}
What this means, is that 3 times the number of stations assigned to district one plus one times the number of stations assigned to district two, plus 2 times the number of stations assigned to district three, must be smaller than 9 stations in total (millions of dollars in the problem was a typo).
We see that part(a) does not meet this requirement. In part(a), we found \(u\left ( 0\right ) =2,u\left ( 1\right ) =0,u\left ( 3\right ) =3\), which means \begin{align*} 3u\left ( 0\right ) +u\left ( 1\right ) +3u\left ( 2\right ) & =3\left ( 2\right ) +0+2\left ( 3\right ) \\ & =12 \end{align*}
Which is larger than \(9\). We need to find again \(u\left ( 0\right ) ,u\left ( 1\right ) ,u\left ( 2\right ) \) which still satisfies part (a) requirements of no more than 3 stations for a district and no more than total of 5 stations, but now with the additional constraint of (1) in place at the same time.
The search was repeated using brute force to first find the combinations that meet this criteria, and then the one with the minimum expected damage was selected.
| district 1 | district 2 | district 3 | cost in millions | \(3u_0+u_1+2u_2\) |
| 0 | 0 | 0 | 4.0 | 0 |
| 0 | 0 | 1 | 3.5 | 2 |
| 0 | 0 | 2 | 3.2 | 4 |
| 0 | 0 | 3 | 2.8 | 6 |
| 0 | 1 | 0 | 3.8 | 1 |
| 0 | 1 | 1 | 3.3 | 3 |
| 0 | 1 | 2 | 3.0 | 5 |
| 0 | 1 | 3 | 2.6 | 7 |
| 0 | 2 | 0 | 3.7 | 2 |
| 0 | 2 | 1 | 3.2 | 4 |
| 0 | 2 | 2 | 2.9 | 6 |
| 0 | 2 | 3 | 2.5 | 8 |
| 0 | 3 | 0 | 3.6 | 3 |
| 0 | 3 | 1 | 3.1 | 5 |
| 0 | 3 | 2 | 2.8 | 7 |
| 1 | 0 | 0 | 2.9 | 3 |
| 1 | 0 | 1 | 2.4 | 5 |
| 1 | 0 | 2 | 2.1 | 7 |
| 1 | 0 | 3 | 1.7 | 9 |
| 1 | 1 | 0 | 2.7 | 4 |
| 1 | 1 | 1 | 2.2 | 6 |
| 1 | 1 | 2 | 1.9 | 8 |
| 1 | 2 | 0 | 2.6 | 5 |
| 1 | 2 | 1 | 2.1 | 7 |
| 1 | 2 | 2 | 1.8 | 9 |
| 1 | 3 | 0 | 2.5 | 6 |
| 1 | 3 | 1 | 2.0 | 8 |
| 2 | 0 | 0 | 2.3 | 6 |
| 2 | 0 | 1 | 1.8 | 8 |
| 2 | 1 | 0 | 2.1 | 7 |
| 2 | 1 | 1 | \(1.6^{\ast }\) | 9 |
| 2 | 2 | 0 | 2.0 | 8 |
| 2 | 3 | 0 | 1.9 | 9 |
| 3 | 0 | 0 | 2.2 | 9 |
We see that the combination with the minimum expected damage is when
The following diagram illustrates the branch and bound graph with the new optimal path now highlighted in the think line.
The code used to part(b) is in the appendix.
code to generate the first table for part(a)
code to generate the first table for part(b)
The state equation is (using indices as subscripts from now on, in all that follows as it is easier to read. Therefore \(x\left ( N\right ) \) is written as \(x_{N}\) and similarly for \(u\left ( N\right ) \))\begin{equation} x_{k+1}=x_{k}-u_{k}\tag{1} \end{equation} The cost function to minimize is\begin{equation} J=\sum _{k=0}^{N-1}\left ( x_{k}-u_{k}\right ) ^{2}+u_{k}^{2}\tag{2} \end{equation}
Now we apply Bellman dynamic equation. This diagram shows the overall process.
The branch cost from \(x_{N-1}\) with one step to go is\[ I\left ( x_{N-1},1\right ) =\min _{u_{N-1}\in \Omega _{N-1}}\left \{ J\left ( x_{N-1},u_{N-1}\right ) +\Psi \left ( x_{N}\right ) \right \} \] \(\Psi \left ( x_{N}\right ) \) is the terminal cost. Removing the terminal cost from the rest of the computation to simplify things and replacing \(J\) in the above from (2) gives\begin{align} I\left ( x_{N-1},1\right ) & =\min _{u_{N-1}\in \Omega _{N-1}}\left \{ J\left ( x_{N-1},u_{N-1}\right ) \right \} \nonumber \\ & =\min _{u_{N-1}\in \Omega _{N-1}}\left \{ \left ( x_{N-1}-u_{N-1}\right ) ^{2}+u_{N-1}^{2}\right \} \tag{3} \end{align}
Taking derivative in order to find optimal \(u_{N-1}^{\ast }\) results in\begin{align*} \frac{dI \left (x_{N-1},1\right )}{u_{N-1}} & =0\\ 2\left ( x_{N-1}-u_{N-1}\right ) \left ( -1\right ) +2u_{N-1} & =0\\ 4u_{N-1}-2x_{N-1} & =0\\ u_{N-1}^{\ast } & =\frac{1}{2}x_{N-1} \end{align*}
Using the above \(u_{N-1}^{\ast }\) in (3) gives\begin{align*} I^{\ast }\left ( x_{N-1},1\right ) & =\left ( x_{N-1}-\frac{1}{2}x_{N-1}\right ) ^{2}+\left ( \frac{1}{2}x_{N-1}\right ) ^{2}\\ & =\frac{1}{2}x_{N-1}^{2} \end{align*}
Going back one step
\begin{align*} I\left ( x_{N-2},2\right ) & =\min _{u_{N-2}\in \Omega _{N-2}}\left \{ J\left ( x_{N-2},u_{N-2}\right ) +I^{\ast }\left ( x_{N-1},1\right ) \right \} \\ & =\min _{u_{N-2}\in \Omega _{N-2}}\left \{ \left ( x_{N-2}-u_{N-2}\right ) ^{2}+u_{N-2}^{2}+\frac{1}{2}x_{N-1}^{2}\right \} \end{align*}
Before taking derivatives, we have to make all the stages to be at \(N-2\), and for this we use the state equation to write \(x_{N-1}\) in terms of \(x_{N-2}\) and the above becomes
\begin{align} I\left ( x_{N-2},2\right ) & =\min _{u_{N-2}\in \Omega _{N-2}}\left \{ \left ( x_{N-2}-u_{N-2}\right ) ^{2}+u_{N-2}^{2}+\frac{1}{2}\left ( x_{N-2}-u_{N-2}\right ) ^{2}\right \} \nonumber \\ & =\min _{u_{N-2}\in \Omega _{N-2}}\left \{ \frac{5}{2}u_{N-2}^{2}-3u_{N-2}x_{N-2}+\frac{3}{2}x_{N-2}^{2}\right \} \tag{4} \end{align}
Now we take derivative to find optimal \(u_{N-2}^{\ast }\) \begin{align*} \frac{dI\left ( x_{N-2},1\right ) }{u_{N-2}} & =0\\ 5u_{N-2}-3x_{N-2} & =0\\ u_{N-2}^{\ast } & =\frac{3}{5}x_{N-2} \end{align*}
We go back to (4) and update with the optimal control found to obtain
\begin{align*} I^{\ast }\left ( x_{N-2},2\right ) & =\frac{5}{2}\left ( \frac{3}{5}x_{N-2}\right ) ^{2}-3\left ( \frac{3}{5}x_{N-2}\right ) x_{N-2}+\frac{3}{2}x_{N-2}^{2}\\ & =\frac{3}{5}x_{N-2}^{2} \end{align*}
Going back one more step
\begin{align*} I\left ( x_{N-3},3\right ) & =\min _{u_{N-3}\in \Omega _{N-3}}\left \{ J\left ( x_{N-3},u_{N-3}\right ) +I^{\ast }\left ( x_{N-2},2\right ) \right \} \\ & =\min _{u_{N-3}\in \Omega _{N-3}}\left \{ \left ( x_{N-3}-u_{N-3}\right ) ^{2}+u_{N-3}^{2}+\frac{3}{5}x_{N-2}^{2}\right \} \end{align*}
Before taking derivatives, we have to make all the states to be at \(N-3\), and for this we use the state equation to write \(x_{N-2}\) in terms of \(x_{N-3}\) and the above becomes
\begin{align} I\left ( x_{N-3},3\right ) & =\min _{u_{N-3}\in \Omega _{N-3}}\left \{ \left ( x_{N-3}-u_{N-3}\right ) ^{2}+u_{N-3}^{2}+\frac{3}{5}\left ( x_{N-3}-u_{N-3}\right ) ^{2}\right \} \nonumber \\ & =\min _{u_{N-2}\in \Omega _{N-2}}\left \{ \frac{13}{5}u_{N-3}^{2}-\frac{16}{5}u_{N-3}x_{N-3}+\frac{8}{5}x_{N-3}^{2}\right \} \tag{5} \end{align}
Now we take derivative to find optimal \(u_{N-3}^{\ast }\) \begin{align*} \frac{dI\left ( x_{N-3},3\right ) }{u_{N-3}} & =0\\ \frac{26}{5}u_{N-3}-\frac{16}{5}x_{N-3} & =0\\ u_{N-3}^{\ast } & =\frac{8}{13}x_{N-3} \end{align*}
We go back to (5) and update the cost to obtain
\begin{align*} I^{\ast }\left ( x_{N-3},3\right ) & =\frac{13}{5}\left ( \frac{8}{13}x_{N-3}\right ) ^{2}-\frac{16}{5}\left ( \frac{8}{13}x_{N-3}\right ) x_{N-3}+\frac{8}{5}x_{N-3}^{2}\\ & =\frac{8}{13}x_{N-3}^{2} \end{align*}
Let us do one more step backward,
\begin{align*} I\left ( x_{N-4},4\right ) & =\min _{u_{N-4}\in \Omega _{N-4}}\left \{ J\left ( x_{N-4},u_{N-4}\right ) +I^{\ast }\left ( x_{N-3},3\right ) \right \} \\ & =\min _{u_{N-4}\in \Omega _{N-4}}\left \{ \left ( x_{N-4}-u_{N-4}\right ) ^{2}+u_{N-4}^{2}+\frac{8}{13}x_{N-3}^{2}\right \} \end{align*}
Before taking derivatives, we have to make all the states to be at \(N-4\), and for this we use the state equation to write \(x_{N-3}\) in terms of \(x_{N-4}\) and the above becomes
\begin{align} I\left ( x_{N-4},4\right ) & =\min _{u_{N-4}\in \Omega _{N-4}}\left \{ \left ( x_{N-4}-u_{N-4}\right ) ^{2}+u_{N-4}^{2}+\frac{8}{13}\left ( x_{N-4}-u_{N-4}\right ) ^{2}\right \} \nonumber \\ & =\min _{u_{N-4}\in \Omega _{N-4}}\left \{ \frac{34}{13}u_{N-4}^{2}-\frac{42}{13}u_{N-4}x_{N-4}+\frac{21}{13}x_{N-4}^{2}\right \} \tag{6} \end{align}
Now we take derivative to find optimal \(u_{N-4}^{\ast }\) \begin{align*} \frac{dI\left ( x_{N-4},4\right ) }{u_{N-4}} & =0\\ \frac{68}{13}u_{N-4}-\frac{42}{13}x_{N-4} & =0\\ u_{N-4}^{\ast } & =\frac{21}{34}x_{N-4} \end{align*}
We go back to (6) and update to obtain
\begin{align*} I^{\ast }\left ( x_{N-4},4\right ) & =\frac{34}{13}\left ( \frac{21}{34}x_{N-4}\right ) ^{2}-\frac{42}{13}\left ( \frac{21}{34}x_{N-4}\right ) x_{N-4}+\frac{21}{13}x_{N-4}^{2}\\ & =\frac{21}{34}x_{N-4}^{2} \end{align*}
This is so much fun, so let us do one more step backward,
\begin{align*} I\left ( x_{N-5},5\right ) & =\min _{u_{N-5}\in \Omega _{N-5}}\left \{ J\left ( x_{N-5},u_{N-5}\right ) +I^{\ast }\left ( x_{N-4},4\right ) \right \} \\ & =\min _{u_{N-5}\in \Omega _{N-5}}\left \{ \left ( x_{N-5}-u_{N-5}\right ) ^{2}+u_{N-5}^{2}+\frac{21}{34}x_{N-4}^{2}\right \} \end{align*}
Before taking derivatives, we have to make all the states to be at \(N-5\), and for this we use the state equation to write \(x_{N-4}\) in terms of \(x_{N-5}\) and the above becomes
\begin{align} I\left ( x_{N-5},5\right ) & =\min _{u_{N-5}\in \Omega _{N-5}}\left \{ \left ( x_{N-5}-u_{N-5}\right ) ^{2}+u_{N-5}^{2}+\frac{21}{34}\left ( x_{N-5}-u_{N-5}\right ) ^{2}\right \} \nonumber \\ & =\min _{u_{N-5}\in \Omega _{N-5}}\left \{ \frac{89}{34}u_{N-5}^{2}-\frac{55}{17}u_{N-5}x_{N-5}+\frac{55}{34}x_{N-5}^{2}\right \} \tag{7} \end{align}
Now we take derivative to find optimal \(u_{N-5}^{\ast }\) \begin{align*} \frac{dI\left ( x_{N-5},5\right ) }{u_{N-5}} & =0\\ \frac{89}{17}u_{N-5}-\frac{55}{17}x_{N-5} & =0\\ u_{N-5}^{\ast } & =\frac{55}{89}x_{N-5} \end{align*}
We go back to (7) and update to obtain
\begin{align*} I^{\ast }\left ( x_{N-5},5\right ) & =\frac{89}{34}\left ( \frac{55}{89}x_{N-5}\right ) ^{2}-\frac{55}{17}\left ( \frac{55}{89}x_{N-5}\right ) x_{N-5}+\frac{55}{34}x_{N-5}^{2}\\ & =\frac{55}{89}x_{N-5}^{2} \end{align*}
Summary table of finding
| \(k\) | \(u^{\ast }\left ( N-k\right ) \) | \(I^{\ast }\left ( x_{N-k},k\right ) \) |
| \(1\) | \(\frac{1}{2}x_{N-1}\) | \(\frac{1}{2}x_{N-1}^{2}\) |
| \(2\) | \(\frac{3}{5}x_{N-2}\) | \(\frac{3}{5}x_{N-2}^{2}\) |
| \(3\) | \(\frac{8}{13}x_{N-3}\) | \(\frac{8}{13}x_{N-3}^{2}\) |
| \(4\) | \(\frac{21}{34}x_{N-4}^{2}\) | \(\frac{21}{34}x_{N-4}^{2}\) |
| \(5\) | \(\frac{55}{89}x_{N-5}\) | \(\frac{55}{89}x_{N-5}^{2}\) |
This is generated using bisection of the Fibonacci sequence, let \(\gamma \left ( k\right ) =\frac{\beta \left ( k\right ) }{\alpha \left ( k\right ) }\) where3 \begin{align*} \beta \left ( k\right ) & =3\beta \left ( k-1\right ) -\beta \left ( k-2\right ) \qquad \\ \beta \left ( 0\right ) & =0\\ \beta \left ( 1\right ) & =1 \end{align*}
And4
\begin{align*} \alpha \left ( k\right ) & =3\alpha \left ( k-1\right ) -\alpha \left ( k-2\right ) \qquad \\ \alpha \left ( 0\right ) & =1\\ \alpha \left ( 1\right ) & =1 \end{align*}
Here is program which generates up to \(k=20\)
which gives
\(\frac{1}{2},\frac{3}{5},\frac{8}{13},\frac{21}{34},\frac{55}{89}, \frac{144}{233},\frac{377}{610},\frac{987}{1597},\frac{2584}{4181}, \frac{6765}{10946},\frac{17711}{28657},\frac{46368}{75025}, \frac{121393}{196418},\frac{317811}{514229},\frac{832040}{1346269}, \frac{2178309}{3524578},\frac{5702887}{9227465}, \frac{14930352}{24157817},\frac{39088169}{63245986},\frac{102334155}{165580141} \)
or numerically
\(\{0.5,0.6,0.6153846153846154,0.6176470588235294,0.6179775280898876, 0.6180257510729614, 0.6180327868852459,0.6180338134001252,0.6180339631667066, 0.618033985017358,0.6180339882053251, 0.6180339886704432,0.618033988738303,0.6180339887482036\}\)
The golden ratio is
\[ \varphi =\frac{1+\sqrt{5}}{2}=1.6180339887482036 \]
Therefore in the limit, for large \(k\)
\[ \fbox{$u^\ast \left ( N-k\right ) =\frac{1}{\varphi }x\left ( N-k\right ) $}\]
Given
\begin{align*} x_{1}\left ( k+1\right ) & =\left [ 1+u^{2}\left ( k\right ) \right ] x_{1}\left ( k\right ) \\ x_{2}\left ( k+1\right ) & =\frac{e^{-u\left ( k\right ) x_{1}\left ( k\right ) }}{x_{1}\left ( k\right ) }+2x_{2}\left ( k\right ) \qquad k=0,1,\cdots ,N \end{align*}
And the goal is to minimize the objective function at the terminal stage \(J=x_{1}\left ( N\right ) +x_{2}\left ( N\right ) \). At stage \(N-1\) with one step to go\begin{equation} I\left ( x\left ( N-1\right ) ,1\right ) =\min _{u\left ( N-1\right ) \in \Omega _{N-1}}\left \{ \Psi \left ( x\left ( N\right ) \right ) \right \} \tag{1} \end{equation} Where \(\Psi \left ( x\left ( N\right ) \right ) =I\left ( x\left ( N\right ) ,0\right ) \). Hence \begin{align*} I\left ( x\left ( N-1\right ) ,1\right ) & =\min _{u\left ( N-1\right ) \in \Omega _{N-1}}\left \{ x_{1}\left ( N\right ) +x_{2}\left ( N\right ) \right \} \\ & =\min _{u\left ( N-1\right ) \in \Omega _{N-1}}\left [ 1+u^{2}\left ( N-1\right ) \right ] x_{1}\left ( N-1\right ) +\frac{e^{-u\left ( N-1\right ) x_{1}\left ( N-1\right ) }}{x_{1}\left ( N-1\right ) }+2x_{2}\left ( N-1\right ) \end{align*}
Therefore \(\frac{dI\left ( x\left ( N-1\right ) \right ) }{du\left ( N-1\right ) }=0\) gives\begin{align*} 0 & =2u\left ( N-1\right ) x_{1}\left ( N-1\right ) -e^{-u\left ( N-1\right ) x_{1}\left ( N-1\right ) }\\ e^{-u\left ( N-1\right ) x_{1}\left ( N-1\right ) } & =2u\left ( N-1\right ) x_{1}\left ( N-1\right ) \end{align*}
This is of the form \(e^{-zx}=2zx\) where \(z\rightarrow u\left ( N-1\right ) \) and \(x\rightarrow x_{1}\left ( N-1\right ) \), which has root at \(z=\frac{0.35173}{x}\) (using Mathematica), hence the control law is\[ \fbox{$u^\ast \left ( N-1\right ) =\frac{0.35173}{x_1\left ( N-1\right ) }$}\]
The following diagram shows the layout of the stages we need to use. There are three stages. \(k=2\) is the terminal stage, and \(k=0\) is the initial stage.
We have
\begin{align*} J & =\min _{k=1,2}x_{2}\left ( k\right ) \\ x_{1}\left ( k+1\right ) & =\min \left \{ x_{1}\left ( k\right ) ,x_{2}\left ( k\right ) \right \} +u\left ( k\right ) \\ x_{2}\left ( k+1\right ) & =x_{1}\left ( k\right ) u\left ( k\right ) \\ x_{1}\left ( 0\right ) & =1\\ x_{2}\left ( 0\right ) & =-1 \end{align*}
One step to go, where \(N=2\)
\begin{align*} I\left ( x\left ( N-1\right ) ,1\right ) & =I\left ( x\left ( 2-1\right ) ,1\right ) \\ & =I\left ( x\left ( 1\right ) ,1\right ) \\ & =\max _{u\left ( 1\right ) \in \Omega _{1}}\left \{ x_{2}\left ( 2\right ) \right \} \end{align*}
But \(x_{2}\left ( 2\right ) =x_{1}\left ( 1\right ) u\left ( 1\right ) \) from the state equation, hence
\[ I\left ( x\left ( 1\right ) ,1\right ) =\max _{u\left ( 1\right ) \in \Omega _{1}}\left \{ x_{1}\left ( 1\right ) u\left ( 1\right ) \right \} \]
We need to find \(u\left ( 1\right ) \) which maximizes \(x_{1}\left ( 1\right ) u\left ( 1\right ) \). Since \(u\left ( k\right ) \in \Omega _{k}=\left [ -M,M\right ] \) then
\[ \fbox{$u^\ast \left ( 1\right ) =M \text{sign}\left ( x\left ( 1\right ) \right ) $}\]
Therefore
\begin{align*} I^{\ast }\left ( x\left ( 1\right ) ,1\right ) & =x_{1}\left ( 1\right ) u^{\ast }\left ( 1\right ) \\ & =x_{1}\left ( 1\right ) M\text{ sign}\left ( x\left ( 1\right ) \right ) \\ & =M\left \vert x_{1}\left ( 1\right ) \right \vert \end{align*}
Now we go one more step backward
\[ I\left ( x\left ( 0\right ) ,2\right ) =\max _{u\left ( 0\right ) \in \Omega _{0}}\min \left \{ x_{2}\left ( 1\right ) ,I^{\ast }\left ( x\left ( 1\right ) ,1\right ) \right \} \]
But from state equation, \(x_{2}\left ( 1\right ) =x_{1}\left ( 0\right ) u\left ( 0\right ) \) and since \(x_{1}\left ( 0\right ) =1\) then \(x_{2}\left ( 1\right ) =u\left ( 0\right ) \) and the above becomes
\begin{equation} I\left ( x\left ( 0\right ) ,2\right ) =\max _{u\left ( 0\right ) \in \Omega _{0}}\min \left \{ u\left ( 0\right ) ,M\left \vert x_{1}\left ( 1\right ) \right \vert \right \} \tag{1} \end{equation}
But \begin{align*} x_{1}\left ( 1\right ) & =\min \left \{ x_{1}\left ( 0\right ) ,x_{2}\left ( 0\right ) \right \} +u\left ( 0\right ) \\ & =\min \left \{ 1,-1\right \} +u\left ( 0\right ) \\ & =-1+u\left ( 0\right ) \end{align*}
Therefore (1) becomes
\[ I\left ( x\left ( 0\right ) ,2\right ) =\max _{u\left ( 0\right ) \in \Omega _{0}}\min \left \{ u\left ( 0\right ) ,M\left \vert -1+u\left ( 0\right ) \right \vert \right \} \]
Let \[ F\left ( u\left ( 0\right ) \right ) =\min \left \{ u\left ( 0\right ) ,M\left \vert -1+u\left ( 0\right ) \right \vert \right \} \]
We need to find \(\max _{u\left ( 0\right ) }\min \left ( F\left ( u\left ( 0\right ) \right ) \right ) \). By making small program and adjusting \(M\), to see all the regions, the following result was found for \(u^{\ast }\left ( 0\right ) \)
| \(M\) range | optimal |
| \(0\leq M\leq \sqrt{2}\) | \(u_{0}^{\ast }=\frac{M}{1+M}\) |
| \(\sqrt{2}<M\) | \(u_{0}^{\ast }=M\) |
The following is a plot of the small program written showing \(F\left ( u\left ( 0\right ) \right ) \) for first case \(0\leq M\leq \sqrt{2}\) and another plot showing the case for \(\sqrt{2}<M\). No other cases were found. The program allows one to move a slider to adjust \(M\) and it finds the maximizing \(u^{\ast }\) for \(F\left ( u\right ) \) at each slider change.
Source code for the above
\begin{align*} \begin{pmatrix} x_{1}\left ( k+1\right ) \\ x_{2}\left ( k+1\right ) \end{pmatrix} & =\overset{A}{\overbrace{\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} }}\begin{pmatrix} x_{1}\left ( k\right ) \\ x_{2}\left ( k\right ) \end{pmatrix} +\overset{B}{\overbrace{\begin{pmatrix} 0\\ 1 \end{pmatrix} }}u\left ( k\right ) \\ J & =\sum _{k=0}^{\infty }2x_{1}^{2}\left ( k\right ) +2x_{1}\left ( k\right ) x_{2}\left ( k\right ) +x_{2}^{2}\left ( k\right ) +3u^{2}\left ( k\right ) \end{align*}
Since \(J\) has the form \(x^{T}Qx+u^{T}Ru\), then we need to find \(Q\) first. (\(Q\) is symmetric), therefore\begin{align*} 2x_{1}^{2}\left ( k\right ) +2x_{1}\left ( k\right ) x_{2}\left ( k\right ) +x_{2}^{2}\left ( k\right ) & =\begin{pmatrix} x_{1}\left ( k\right ) & x_{2}\left ( k\right ) \end{pmatrix}\begin{pmatrix} q_{11} & q_{12}\\ q_{12} & q_{22}\end{pmatrix}\begin{pmatrix} x_{1}\left ( k\right ) \\ x_{2}\left ( k\right ) \end{pmatrix} \\ & =\begin{pmatrix} x_{1}\left ( k\right ) q_{11}+x_{2}\left ( k\right ) q_{12} & x_{1}\left ( k\right ) q_{12}+x_{2}\left ( k\right ) q_{22}\end{pmatrix}\begin{pmatrix} x_{1}\left ( k\right ) \\ x_{2}\left ( k\right ) \end{pmatrix} \\ & =x_{1}^{2}\left ( k\right ) q_{11}+2q_{12}x_{1}\left ( k\right ) x_{2}\left ( k\right ) +x_{2}^{2}\left ( k\right ) q_{22} \end{align*}
Comparing coefficients, we see that \(q_{11}=2,q_{22}=1,q_{12}=1\), hence\[ Q=\begin{pmatrix} 2 & 1\\ 1 & 1 \end{pmatrix} \] And \(R\) is scalar\[ R=3 \] Therefore, the discrete algebraic Riccati equation is\begin{equation} P=APA^{T}-A^{T}PB\left ( R+B^{T}PB\right ) ^{-1}B^{T}PA+Q \tag{1} \end{equation} Small note: The above is what we had in lecture notes. Matlab has the above in its help pages slightly different. The first term is written as \(A^{T}PA\) instead of \(APA^{T}\) as we had.
Using Matlab
\[ P=\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix} \] Notice that \(P\) is symmetric as expected.
Let us check it is correct first. Substituting \(P\) in RHS of (1) gives\begin{align*} P & =\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix}\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix}\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} ^{T}-\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} ^{T}\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix}\begin{pmatrix} 0\\ 1 \end{pmatrix} \\ & \left ( 3+\begin{pmatrix} 0\\ 1 \end{pmatrix} ^{T}\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix}\begin{pmatrix} 0\\ 1 \end{pmatrix} \right ) ^{-1}\begin{pmatrix} 0\\ 1 \end{pmatrix} ^{T}\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix}\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix} \\ & +\begin{pmatrix} 2 & 1\\ 1 & 1 \end{pmatrix} \end{align*}
The above gives \(\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix} \) which is \(P\). Yes, it satisfies DARE. Now, using\begin{align*} u^{\ast } & =-\left ( R+B^{T}PB\right ) ^{-1}B^{T}PAx\\ & =-\left ( 3+\begin{pmatrix} 0\\ 1 \end{pmatrix} ^{T}\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix}\begin{pmatrix} 0\\ 1 \end{pmatrix} \right ) ^{-1}\begin{pmatrix} 0\\ 1 \end{pmatrix} ^{T}\begin{pmatrix} 3.7841 & 1.6815\\ 1.6815 & 4.4022 \end{pmatrix}\begin{pmatrix} 0 & 1\\ 1 & 0 \end{pmatrix}\begin{pmatrix} x_{1}\left ( k\right ) \\ x_{2}\left ( k\right ) \end{pmatrix} \\ & =-0.59472x_{1}\left ( k\right ) -0.22716x_{2}\left ( k\right ) \end{align*}
Therefore, comparing the above to \(u^{\ast }\left ( k\right ) =K_{1}x_{1}\left ( k\right ) +K_{2}x_{2}\left ( k\right ) \) we see that \begin{align*} K_{1} & =-0.59472\\ K_{2} & =-0.22716 \end{align*}
To simulate the result under the new control law \(u^{\ast }\), we need some initial condition on state. Below is small script which simulate this up to \(k=20\) stages with the plot generated
