Overview

This project contains supplemental Matlab code for the article:

M. Thelander Andrén, B. Bernhardsson, A. Cervin and K. Soltesz, "On Event-Based Sampling for LQG-Optimal Control", In Proc. 56th IEEE Conf. on Decision and Control, 2017

It demonstrates a numerical method for computing the optimal event-based sampling scheme for the continious-time LQG problem. The problem is related to an elliptic, convection-diffusion type of partial-differential equation (PDE) with free boundary, a so called Stefan problem. The PDE is:

	\forall x_{\text{\tiny H}}\in \mathbb{R}^n: \begin{cases}
	x_{\text{\tiny H}}^\intercal Qx_{\text{\tiny H}} - J + x_{\text{\tiny H}}^\intercal A^\intercal \nabla V + \frac{1}{2}\text{Tr}(R\nabla^2V) = \frac{\partial V}{\partial t}, \\  
	V(x_{\text{\tiny H}},t)\leq 0, 
	\end{cases}\quad\quad
	 \forall x_{\text{\tiny H}} \in \partial \Omega: 
	\begin{cases}
	V(x_{\text{\tiny H}},t) = 0,\\
	\nabla V = 0.
	\end{cases}

The solution to this PDE is the value function V. In stationarity, the free boundary \partial \Omega is the threshold on the vector state x_{\text{\tiny H}} which defines the optimal sampling scheme. The current implementation handles second-order systems. For more details, we refer to our article.

Running the Code

For getting started, we recommend running and modifying the script script_optimal_bound_2D.m. The stationary solution of V is computed by running the function optimal_bound_2D, found in functions/optimal_bound_2D.m. A typical run is:

%Define system
A = [-1 20; 20  -1];
q11 = 1;
q22 = 1;
J=1;

%Gridding
dx = 3e-2;
x1min=-2;
x1max = 2;
x2min= -2;
x2max = 2;
x1vec = x1min:dx:x1max;
x2vec = x2min:dx:x2max;
h = 5e-4;

%Enable/disable plotting and printing
plot_progress = 0;
print_progress = 1;

%Compute V
V = optimal_bound_2d(A, q11, q22, J, x1vec, x2vec, h, plot_progress, print_progress);

Without loss of generality, optimal_bound_2D assumes that the system is of a form where R = I and Q = \text{diag}(q_{11}, q_{22}) (the system state can always be transformed to achieve this). The input arguments of optimal_bound_2D are:

• A : The 2\times2 system matrix of the linear reset system.
• q11,q22 : The diagonal elements of the (diagonal) cost matrix Q.
• J : The total cost of state deviation and sampling.
• x1vec,x2vec : The grid vectors of the two state dimensions. Each vector is assumed to have equidistantly spaced elements.
• h : The time step of the simulation.
• plot_progress : Specifies if V and \partial \Omega should be plotted during the simulation (plot_progress = 1) or not (plot_progress = 0).
• print_progress : Specifies if iteration info should be printed to the command window during the simulation (print_progress = 1) or not (print_progress = 0).

The output argument V is a matrix with dimensions compatible with x1vecand x2vec. For instance, you can call surf(x1vec, x2vec, V) to generate a surface plot of V, or contour(x1vec, x2vec, V, [0 0]) to plot the bound \partial \Omega.

Plotting during Simulation

With the input argument plot_progress = 1, a surface plot of V and a plot of \partial \Omega will be updated every 10th iteration of the simulation.

Printing during Simulation

With the input argument print_progress = 1, information of every 10th iteration will be printed to the command window during the simulation. An example is shown below:

Iterations : 10 out of 50000 | Maximum Vdiff = 0.00049488 | V(0) = -0.0049714

The printout displays the current iteration number (in this case 10) out of the maximum number of iterations (50000). The Maximum Vdiff value is the largest absolute value between the current and last iteration V, and is the basis for determining when stationarity is reached. The value V(0) is the current value of V at the origin. It is linked to the per-sample cost \rho, as V(0) = -\rho (see the article).

About the Numerical Method

The numerical solution method is based on simulating the non-stationary version of the PDE using finite-difference approximations of the differential operators and gridding in time and space. In this implementation, we use the backward-time central-space (BTCS) finite-difference approximation (see e.g. this article for more details).

We simulate the PDE by iterating the BTCS scheme, while at each iteration enforcing the inequality V \leq 0 by the assignment V = min(V,0). This has the effect of moving the free boundary \partial \Omega over time. The simulation is then progressed until a stopping criterion based on stationarity of V is met, or if the maximum number of iterations is reached. The stopping criterion is an upper threshold on the maximum absolute difference of V between to successive iterations.

The current implementation is for second-order systems only, i.e. x_{\text{\tiny H}} \in \mathbb{R}^2. However, this method of applying BTCS can be directly generalized to higher order systems as well, although at a higher computational cost (the time complexity scales exponentially with the state dimension).