{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Part1: Constrained Least Squares Problems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

1. Linearly constrained least squares

\n", "

\n", "\n", "the (linearly) constrained least squares problem (CLS) is\n", "\n", "minimize $\\|A x-b\\|^{2}$\n", "\n", "subject to $C x=d$\n", "\n", "$\\|A x-b\\|^{2}$ is the objective function and Cx=d are equality constraints, x is an n-vector, A is a m × n matrix, b is an m-vector, C is a p × n matrix, and d is a p-vector.\n", "\n", "$\\hat{\\mathbf{x}}$ is a solution of CLS if $C \\hat{x}=d$ and $\\|A \\hat{x}-b\\| \\leq\\|A x-b\\|^{2}$ holds for any n-vector x that satisfies $C x=d$\n", "\n", "\n", "

2.Methods of solving the constrained least squares problem

\n", "

\n", "2.1 KKT equations\n", "optimality conditions in matrix-vector form:\n", "\n", "$$2\\left(A^{T} A\\right) \\hat{x}-2 A^{T} \\mathbf{b}+C^{T} \\mathbf{z}=0, C \\hat{\\mathbf{x}}=\\mathbf{d}$$\n", "put these together to get KKT conditions $\\left[\\begin{array}{cc}{2\\left(A^{T} A\\right)} & {C^{T}} \\\\ {C} & {0}\\end{array}\\right]\\left[\\begin{array}{l}{\\hat{x}} \\\\ {\\hat{z}}\\end{array}\\right]=\\left[\\begin{array}{c}{2 A^{T} b} \\\\ {d}\\end{array}\\right]$\n", "\n", "then we get \n", "$\\left[\\begin{array}{c}{\\hat{\\mathrm{x}}} \\\\ {\\hat{\\mathrm{z}}}\\end{array}\\right]=\\left[\\begin{array}{cc}{2\\left(A^{T} A\\right)} & {C^{T}} \\\\ {C} & {0}\\end{array}\\right]^{-1}\\left[\\begin{array}{c}{2 A^{T} \\mathrm{b}} \\\\ {\\mathrm{d}}\\end{array}\\right]$\n", "\n", "2.2 QR factorization\n", "\n", "

3.Least norm problem (an important special case of Least Squares Problems)

\n", "

\n", "minimize $\\|x\\|^{2}$\n", "\n", "subject to $c_{i}^{T} x=d_{i}, \\mathrm{i}=1, \\ldots, \\mathrm{p}$\n", "\n", "3.1 form Lagrangian function, with Lagrange multipliers $ z_{1}, \\ldots, z_{p} . \\mathrm{L}(\\mathrm{x}, \\mathrm{z})=\\mathrm{f}(\\mathrm{x})+z_{1}\\left(c_{1}^{T} x-d_{1}\\right)+\\cdots+z_{p}\\left(c_{p}^{T} x-d_{p}\\right)$\n", "\n", "3.2 optimality conditions are $\\frac{\\partial \\mathrm{L}}{\\partial x_{i}}(\\hat{x}, \\hat{z})=2 \\sum_{j=1}^{n}\\left(A^{T} A\\right)_{i j} \\hat{x}_{j}-2\\left(A^{T} b\\right)_{i}+\\sum_{j=1}^{p} \\hat{z}_{j}\\left(c_{j}\\right)_{i}=0, i=1, \\ldots, n, \\frac{\\partial \\mathrm{L}}{\\partial z_{i}}(\\hat{x}, \\hat{z})=c_{i}^{T} x-d_{i}=0, \\mathrm{i}=$\n", "$1, \\ldots, \\mathrm{p},$\n", "\n", "3.3 Lagrange Multipliers = Derivatives of the Cost\n", "\n", "Our first example is in two dimensions. The function F is quadratic. The set K is linear.\n", " \n", "$$F(x) = x_{1}^2 + x_{2}^2$$ on the line $$K: a_{1}x_{1} + a_{2}x_{2} = b$$\n", " \n", "On the line K, we are looking for the point that is nearest to (0,0). The cost F(x) is distance squared.\n", "We can discover this from simple calculus, after we bring the constraint equation $a_{1}x_{1} + a_{2}x_{2} = b$ into the function $F(x) = x_{1}^2 + x_{2}^2$\n", "\n", "\n", " Multiply $a_{1}x_{1} + a_{2}x_{2} - b$ by an unknown multiplier λ and add it to F(x)\n", " \n", " Lagrangian $L(x,λ) = F(x) + λ(a_{1}x_{1} + a_{2}x_{2} - b)$\n", "\n", " $= x_{1}^2 + x_{2}^2 + λ(a_{1}x_{1} + a_{2}x_{2} - b)$\n", " \n", " Set the derivatives $∂L/∂x_{1}$ and $∂L/∂x_{2}$ and $ ∂L/∂λ$ to zero.\n", " Solve those three equations for $x_{1},x_{2},λ$\n", "\n", "$$∂L/∂x_{1} = 2x_{1} + λa_{1} = 0$$\n", "$$∂L/∂x_{2} = 2x_{2} + λa_{2} = 0$$\n", "$$∂L/∂λ = a_{1}x_{1} + a_{2}x_{2} - b$$\n", "\n", "\n", "We can get $$λ = -2b/(a_{1}^2 + a_{2}^2)$$\n", "$$x_{1} = -λa_{1}/2 = a_{1}b/(a_{1}^2 + a_{2}^2)$$\n", "$$x_{2} = -λa_{2}/2 = a_{2}b/(a_{1}^2 + a_{2}^2)$$\n", "$$x_{1}^2 + x_{2}^2 = b^2/(a_{1}^2 + a_{2}^2)$$\n", "\n", "The derivative of the minimum cost with respect to the comtrainst level b is minus the Lagrange multiplier:\n", "$$d/db(b^2/(a_{1}^2 + a_{2}^2)) = 2b/(a_{1}^2 + a_{2}^2) = -λ$$\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part2: Solving constrained least problems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

1. Constraint least squares via KKT equations

\n", "

\n", "Algorithm\n", "

\n", "

\n", "Step 1 : Form Gram matrix. Compute $A^{T} A$\n", "

\n", "Step 2 : Solve KKT equations by QR factorization and back substitution.\n", "

\n", "Time complexity : $m n^{2}+2(n+p)^{3}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "

\n", "An example of Constraint least squares via KKT equations in Julia: \n", "

\n", "Firstly, we randomly generate A, b, C, d" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "m = 10; n = 5; p = 2;\n", "A = randn(m,n); b = randn(m); C = randn(p,n); d = randn(p);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Define the function cls_solve_kkt() to implement the algorithm." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "function cls_solve_KKT(A,b,C,d)\n", "m, n = size(A)\n", "p, n = size(C)\n", "G = A'*A # Gram matrix\n", "KKT = [2*G C'; C zeros(p,p)] # KKT matrix\n", "xzhat = KKT \\ [2*A'*b; d]\n", "return xzhat[1:n,:]\n", "end;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test the function cls_solve_kkt()." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5×1 Array{Float64,2}:\n", " 0.0990152352422455 \n", " 0.11870067836607096\n", " -0.490651246965513 \n", " 0.6001981250682451 \n", " -0.22194625122020822" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cls_solve_KKT(A,b,C,d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "

\n", "

\n", "

\n", "

2. Constraint least squares via QR factorization

\n", "

\n", "Algorithm\n", "\n", "Step 1 : Compute the QR factorizations. $\\left[\\begin{array}{c}{A} \\\\ {C}\\end{array}\\right]=\\left[\\begin{array}{l}{Q_{1}} \\\\ {Q_{2}}\\end{array}\\right] R, \\quad Q_{2}^{T}=\\tilde{Q} \\tilde{R}$\n", "

\n", "Step 2 : Compute $\\tilde{R}^{-T} d$ by forward substitution.\n", "

\n", "Step 3 : Form right-hand side and solve $\\tilde{R} w=2 \\tilde{Q}^{T} Q_{1}^{T} b-2 \\tilde{R}^{-T} d$ via back substitution.\n", "

\n", "Step 4 : Compute$\\hat{x}$ Form right-hand side and solve $R \\hat{x}=Q_{1}^{T} b-(1 / 2) Q_{2}^{T} w$ by back substitution.\n", "

\n", "Time complexity : $2(m+p) n^{2}+2 n p^{2}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "

\n", "An example of Constraint least squares via QR factorization in Julia: \n", "

\n", "\n", "Define the function cls_solve_QR() to implement the algorithm." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra\n", "function cls_solve_QR(A,b,C,d)\n", "m, n = size(A)\n", "p, n = size(C)\n", "Q, R = qr([A; C])\n", "Q = Matrix(Q)\n", "Q1 = Q[1:m,:]\n", "Q2 = Q[m+1:m+p,:]\n", "Qtil, Rtil = LinearAlgebra.qr(Q2')\n", "Qtil = Matrix(Qtil)\n", "w = Rtil \\ (2*Qtil'*Q1'*b - 2*(Rtil'\\d))\n", "return xhat = R \\ (Q1'*b - Q2'*w/2)\n", "end;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test the function cls_solve_QR() and compare the result with cls_solve_KKT(). \n", "

\n", "Must be same!" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Array{Float64,1}:\n", " 0.0990152352422457 \n", " 0.11870067836607087\n", " -0.4906512469655131 \n", " 0.6001981250682452 \n", " -0.22194625122020828" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cls_solve_QR(A,b,C,d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "

\n", "

\n", "3. Sparse constrained least squares" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To deal with sparse matrices, the simplest way is to replace the QR factorizations with sparse QR factorizations in the previous algorithms. Fortunately, the built-in function qr() can also realize sparse QR factorizations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "

\n", "An example of solving sparse constraint least squares in Julia: \n", "

\n", "Notice that unlike cls_solve_KKT(), this function assumes b and d are vectors. The following formulation will generate a sparse set of equations to solve if A and C are sparse. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "function cls_solve_sparse(A,b,C,d)\n", "m, n = size(A)\n", "p, n = size(C)\n", "bigA = [ zeros(n,n) A' C';\n", "A -I/2 zeros(m,p) ;\n", "C zeros(p,m) zeros(p,p) ]\n", "xyzhat = bigA \\ [zeros(n) ; b ; d]\n", "return xhat = xyzhat[1:n]\n", "end;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A random formulation of A, b, C, d." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "m = 100; n = 50; p = 10;\n", "A = randn(m,n); b = randn(m); C = randn(p,n); d = randn(p);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, compare based on result of the algorithms of computing via KKT equation and QR factorization. We found the two algorithms agree." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8.863201352429109e-15" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x1 = cls_solve_KKT(A,b,C,d);\n", "x2 = cls_solve_sparse(A,b,C,d);\n", "norm(x1-x2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "

\n", "

\n", "4. Solution of least norm problem\n", "

\n", "For the least norm problem, the KKT equation is reduced to\n", "\\begin{equation}\n", "\\left[\\begin{array}{cc}{2 I} & {C^{T}} \\\\ {C} & {0}\\end{array}\\right]\\left[\\begin{array}{l}{\\hat{x}} \\\\ {\\hat{z}}\\end{array}\\right]=\\left[\\begin{array}{l}{0} \\\\ {d}\\end{array}\\right]\n", "\\end{equation}\n", "\n", "Step 1 : QR factorization. Compute the QR factorization $C^{T}=Q R$\n", "

\n", "Step 2 : Compute$\\hat{x}$. Solve $R^{T} y=d$ by forward substitution.\n", "

\n", "Step 3 : Compute $\\hat{x}=Q y$\n", "

\n", "\n", "Time complexity :$2 n p^{2}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "

\n", "Comparison of sloving least norm problem via different methods in Julia: \n", "

\n", "In julia, the backlash operator can be used to find the solution of equation $C x=d$.\n", "Here is the comparison of the results of solving a least norm problem by different methods." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6.832042562912264e-15" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p = 50; n = 500;\n", "C = randn(p,n); d = randn(p);\n", "# Solve using backslash\n", "x1 = C\\d; \n", "# Solve using cls_solve, which uses KKT system\n", "x2 = cls_solve_KKT(Matrix{Float64}(I, n, n), zeros(n), C, d);\n", "# Using pseudo-inverse\n", "x3 = pinv(C)*d; \n", "norm(x1-x2)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5.904483805221459e-16" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(x1-x3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is obvious that the three methods agree when they deal with least norm problem. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part3: Constrained least squares applications" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two main applications for constrained least squares, one is Linear quadratic control and another is Linear quadratic state estimation. In the control problem, we can choose the inputs; they are under our control. Once we choose the inputs, we know the state sequence. The inputs are typically actions that we take to affect the state trajectory. In the estimation problem, the inputs (called process noise in the estimation problem) are unknown, and the problem is to guess them." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Application 1: Linear quadratic control \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# dynamics equations\n", "$$x_{t+1}=A_{t} x_{t}+B_{t} u_{t}$$\n", "$$y_{t}=C_{t} x_{t}$$\n", "\n", "$$\n", "\\begin{aligned} J_{\\text {output }} &=\\left\\|y_{1}\\right\\|^{2}+\\cdots+\\left\\|y_{T}\\right\\|^{2}=\\left\\|C_{1} x_{1}\\right\\|^{2}+\\cdots+\\left\\|C_{T} x_{T}\\right\\|^{2} \\\\ J_{\\text {input }} &=\\left\\|u_{1}\\right\\|^{2}+\\cdots+\\left\\|u_{T-1}\\right\\|^{2} \\end{aligned}\n", "$$\n", "\n", "The linear quadratic control problem (with initial and final state constraints) is\n", "\n", "minimize $J_{\\text {output }}+\\rho J_{\\text {input }}$ -> $\\|\\tilde{A} z-\\tilde{b}\\|^{2}$\n", "\n", "subject to $x_{t+1}=A_{t} x_{t}+B_{t} u_{t}$\n", "\n", " and $x_{1}=x^{\\text {init }}, \\quad x_{T}=x^{\\text {des }}$\n", " \n", "creating vector z which includes all the variable $z=\\left(x_{1}, \\ldots, x_{T}, u_{1}, \\ldots, u_{T-1}\\right)$\n", " \n", " \n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Example" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "H = randn(2,2); # creating 2*2 matrix" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra\n", "identity_matrix=Matrix{Float64}(I, 3, 3);# creating 3*3 identity matrix" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "6×6 Array{Float64,2}:\n", " 1.04291 0.0107771 0.0 0.0 0.0 0.0 \n", " -0.534264 0.578265 -0.0 0.0 -0.0 0.0 \n", " 0.0 0.0 1.04291 0.0107771 0.0 0.0 \n", " -0.0 0.0 -0.534264 0.578265 -0.0 0.0 \n", " 0.0 0.0 0.0 0.0 1.04291 0.0107771\n", " -0.0 0.0 -0.0 0.0 -0.534264 0.578265 " ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kron(identity_matrix,H) # using kronecker function to block diagonal matrix" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "function cls_solve(A,b,C,d)\n", "m, n = size(A)\n", "p, n = size(C)\n", "Q, R = qr([A; C])\n", "Q = Matrix(Q)\n", "Q1 = Q[1:m,:]\n", "Q2 = Q[m+1:m+p,:]\n", "Qtil, Rtil = qr(Q2')\n", "Qtil = Matrix(Qtil)\n", "w = Rtil \\ (2*Qtil'*Q1'*b - 2*(Rtil'\\d))\n", "return xhat = R \\ (Q1'*b - Q2'*w/2)\n", "end;" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "function eye(k)\n", " matrix=I+zeros(k,k)\n", "return matrix\n", "end;" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "function lqr(A,B,C,x_init,x_des,T,rho)\n", "n = size(A,1)\n", "m = size(B,2)\n", "p = size(C,1)\n", "q = size(x_init,2)\n", "Atil = [ kron(eye(T), C) zeros(p*T,m*(T-1)) ;zeros(m*(T-1), n*T) sqrt(rho)*eye(m*(T-1)) ]\n", "btil = zeros(p*T + m*(T-1), q)\n", "# We’ll construct Ctilde bit by bit\n", "Ctil11 = [ kron(eye(T-1), A) zeros(n*(T-1),n) ] -[ zeros(n*(T-1), n) eye(n*(T-1)) ]\n", "Ctil12 = kron(eye(T-1), B)\n", "Ctil21 = [eye(n) zeros(n,n*(T-1)); zeros(n,n*(T-1)) eye(n)]\n", "Ctil22 = zeros(2*n,m*(T-1))\n", "Ctil = [Ctil11 Ctil12; Ctil21 Ctil22]\n", "dtil = [zeros(n*(T-1), q); x_init; x_des]\n", "z = cls_solve(Atil,btil,Ctil,dtil)\n", "x = [z[(i-1)*n+1:i*n,:] for i=1:T]\n", "u = [z[n*T+(i-1)*m+1 : n*T+i*m, :] for i=1:T-1]\n", "y = [C*xt for xt in x]\n", "return x, u, y\n", "end;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$x=[ x[1],x[2],x[3],x[4],x[5],......,x[T]]$\n", "\n", "$u=[ u[1],u[2],u[3],u[4],u[5],......,u[T]]$\n", "\n", "$y=[ y[1],y[2],y[3],y[4],y[5],......,y[T]]$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "25\n", "\n", "\n", "50\n", "\n", "\n", "75\n", "\n", "\n", "100\n", "\n", "\n", "0.0\n", "\n", "\n", "0.1\n", "\n", "\n", "0.2\n", "\n", "\n", "0.3\n", "\n", "\n", "0.4\n", "\n", "\n", "\n" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using LinearAlgebra\n", "A = [ 0.855 1.161 0.667;\n", "0.015 1.073 0.053;\n", "-0.084 0.059 1.022 ];\n", "B = [-0.076; -0.139; 0.342 ];\n", "C = [ 0.218 -3.597 -1.683 ];\n", "n = 3; p = 1; m = 1;\n", "x_init = [0.496; -0.745; 1.394];\n", "x_des = zeros(n,1);\n", "T = 100;\n", "yol = zeros(T,1);\n", "Xol = [ x_init zeros(n, T-1) ];\n", "for k=1:T-1\n", "Xol[:,k+1] = A*Xol[:,k];\n", "end;\n", "yol = C*Xol;\n", "using Plots\n", "plot(1:T, yol', legend = false)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7738942551160125" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rho = 0.2;\n", "T = 100;\n", "x, u, y = lqr(A,B,C,x_init,x_des,T,rho)\n", "J_input = norm(u)^2" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3.7829986463324596" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "J_output = norm(y)^2" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "25\n", "\n", "\n", "50\n", "\n", "\n", "75\n", "\n", "\n", "100\n", "\n", "\n", "-0.2\n", "\n", "\n", "-0.1\n", "\n", "\n", "0.0\n", "\n", "\n", "0.1\n", "\n", "\n", "0.2\n", "\n", "\n", "0.3\n", "\n", "\n", "t\n", "\n", "\n", "u_t\n", "\n", "\n", "\n" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(1:T-1, vcat(u...), legend = false, xlabel=\"t\",ylabel= \"u_t\")" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "\n", "\n", "25\n", "\n", "\n", "50\n", "\n", "\n", "75\n", "\n", "\n", "100\n", "\n", "\n", "0.0\n", "\n", "\n", "0.1\n", "\n", "\n", "0.2\n", "\n", "\n", "0.3\n", "\n", "\n", "0.4\n", "\n", "\n", "t\n", "\n", "\n", "y_t\n", "\n", "\n", "\n" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "plot(1:T, vcat(y...), legend=false, xlabel = \"t\",ylabel = \"y_t\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Application2: Linear quadratic state estimation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear dynamical system equations\n", "$x_{t+1}=A_{t} x_{t}+B_{t} w_{t}$\n", "\n", "$y_{t}=C_{t} x_{t}+v_{t}$\n", "\n", "$J_{\\mathrm{meas}}=\\left\\|v_{1}\\right\\|^{2}+\\cdots+\\left\\|v_{T}\\right\\|^{2}=\\left\\|C_{1} x_{1}-y_{1}\\right\\|^{2}+\\cdots+\\left\\|C_{T} x_{T}-y_{T}\\right\\|^{2}$\n", "\n", "$J_{\\mathrm{proc}}=\\left\\|w_{1}\\right\\|^{2}+\\cdots+\\left\\|w_{T-1}\\right\\|^{2}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Least squares state estimation. We will make our guesses of x1, . . . , xT and\n", "w1, . . . , wT −1 so as to minimize a weighted sum of our objectives, subject to the dynamics constraints:\n", " \n", " minimize $J_{\\mathrm{meas}}+\\lambda J_{\\mathrm{proc}}$\n", "\n", " subject to $x_{t+1}=A_{t} x_{t}+B_{t} w_{t}$" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.0.3", "language": "julia", "name": "julia-1.0" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.0.3" } }, "nbformat": 4, "nbformat_minor": 2 }