Example_network_estimation

Example 4: Estimation Over a Network

In this example, a physical plant, its output digitally transmitted through a network, and a state estimator are modeled in Simulink as a hybrid system.

Contents

• Mathematical Model
• Steps to Run Model
• Simulink Model
• Example Output
• References

The files for this example are found in the package hybrid.examples.network_estimation :

• initialize.m
• network_est.slx
• postprocess.m

The contents of this package are located in Examples\+hybrid\+examples\network_estimation (clicking this link changes your working directory).

Mathematical Model

Consider a physical process given in terms of the state-space model

$$ \begin{array}{l} \dot{x}=A x,\\\ y= M x \end{array} $$

where $x \in \mathbf{R}^{n_P}$ is the state and $y \in \mathbf{R}^{r_P}$ is the measured output. The output is digitally transmitted through a network. At the other end of the network, a computer receives the information and runs an algorithm that takes the measurements of $y$ to estimate the state $x$ of the physical process.

We consider an estimation algorithm with a state variable $\hat{x} \in \mathbf{R}^{n_P}$ , which denotes the estimate of $x$ , that is appropriately reset to a new value involving the information received. More precisely, denoting the transmission times by $t_i$ and $L$ a constant matrix to be designed, the estimation algorithm updates its state as follows

$$ \hat{x}^+=\hat{x}+L(y-M\hat{x}) $$

at every instant information is received. In between such events, the algorithm updates its state continuously so as to match the evolution of the state of the physical process, that is, via

$$ \dot{\hat{x}} = A \hat{x} $$

Modeling the network as a hybrid system, which, in particular, assumes zero transmission delay, the state variables of the entire system are $j_N$ , $\tau_N \in \mathbf{R}$ , $m_N \in \mathbf{R}^{r_P}$ , and $x,\ \hat{x}\in \mathbf{R}^{n_P}$ . Then, transmissions occur when $\tau_N \leq 0$ , at which events the state of the network is updated via

$$ \tau_N^+ \in [T^{*\min}_{N},T^{*\max}_{N}],\quad j_N^+ = j_N+1, \qquad m_N^+ = y $$

and the state of the algorithm is updated via $\hat{x}^+=\hat{x}+L(y-M\hat{x})$ . Note that since the state of the physical process does not change at such events, we can use the following trivial difference equation to update it at the network events:

$$ x^+ = x $$

In between events, the state of the network is updated as

$$ \dot{\tau}_N = -1, \qquad \dot{j}_N,\dot{m}_N = 0 $$

the state of the algorithm changes continuously according to $\dot{\hat{x}} = A \hat{x}$ , and the state of the physical process changes according to $\dot{x}=A x , y= M x$ . Based on the equations above, we pick the following data for each of the subsystems in the Simulink Model:

Physical process:

$$ \begin{array}{l@{}l} f_P(x, u):= Ax + Bu, \quad &C_P := \mathbf{R}^{n_p}\times\mathbf{R}, \\\ G_P(x, u):= x, \quad &D_P := \emptyset,\\\ y = h(x) : =Mx \end{array}$$

where

$$\begin{array}{l@{}l} A = \left[\begin{tabular}{cccc} 0&1&0&0\\\ -1&0&0&0\\\ -2&1&-1&0\\\ 2&-2&0&-2 \end{tabular}\right],\quad B = \left[\begin{tabular}{c} 0\\0\\1\\0 \end{tabular}\right],\quad M = \left[\begin{tabular}{cccc}1&0&0&0\end{tabular}\right], \end{array}$$

$n_p=4$ , $r_p=1$ , $x \in\mathbf{R}^{n_p}$ , $y\in\mathbf{R}^{r_p}$ , and $u\in\mathbf{R}$ .

Network:

$$\begin{array}{l@{}l} f(x_N, u_N):= \left[\begin{tabular}{c} 0\\0\\-1 \end{tabular}\right],\quad &C := \{(x_N,u_N)\mid \tau_N \geq 0 \},\\\ g(x_N, u_N):= \left[\begin{tabular}{c} $u_N$ \\ $j_N +1$ \\ $\tau_r$ \end{tabular}\right],\quad &D :=\{(x_N,u_N)\mid \tau_N \leq 0 \},\\\ y_N = h(x_N) : =x_N \end{array}$$

where

$\tau_r\in[T^{*\min}_{N},T^{*\max}_{N}]$ is a random variable that models the time in-between communication instances. Also, the sate and the input are given by $x_N=(m_N,j_N,\tau_N)\in\mathbf{R}\times\mathbf{R}\times\mathbf{R}$ , and $u_N=y\in\mathbf{R}^{r_p}$ , respectively.

Estimator:

$$\begin{array}{ll} f(x_E, u_E):= \left[\begin{tabular}{cc} $A$ &0 \\ 0 & 0 \end{tabular}\right]x_E + \left[\begin{tabular}{c} $B$\\0 \end{tabular}\right]v,\quad &C := \{(x_E,u_E) \mid j_E = j_N \},\\\ g((\hat{x},j_E), u_E):= \left[\begin{tabular}{c} $\hat{x} + L(m_N-M\hat{x})$\\\ $j_N$ \end{tabular}\right],\quad &D :=\{(x,u) \mid j_E\in\mathbf{N}\setminus j_N\}, \end{array}$$

where $L$ , which is designed as in [1], is given by

$$\begin{array}{l} L := \left[\begin{tabular}{c} 1 \\ -0.9433\\ -0.6773\\1.6274 \end{tabular}\right], \end{array}$$

the input and the state are given by $u_E = (x_N,v) = ((m_N,j_N,\tau_N),v)\in\mathbf{R}\times\mathbf{R}\times\mathbf{R}\times\mathbf{R}$ , and $x_E=(\hat{x},j_E)\in\mathbf{R}^{4}\times\mathbf{R}$ , respectively. Notice that the estimator block (Estimator) in the Simulink Model is implemented using a regular hybrid system block with embedded functions.

For each hybrid system in the Simulink Model (Continuous Process, network, and Estimator) we have the following Matlab embedded functions that describe the sets $C$ and $D$ and functions $f$ and $g$ .

Steps to Run Model

The following procedure is used to simulate the example using the model in the file network_est.slx :

• Open network_est.slx in Simulink.
• Double-click the block labeled Double Click to Initialize .
• To start the simulation, click the run button or select Simulation>Run .
• Once the simulation finishes, click the block labeled Double Click to Plot Solutions . Several plots of the computed solution will open.

Simulink Model

The following diagram shows the Simulink model of the estimator over a network example. The contents of the blocks flow map f , flow set C , etc., are shown below. When the Simulink model is open, the blocks can be viewed and modified by double clicking on them.

The Simulink blocks for the hybrid systems in this example are included below.

Continuous Process:

flow map f block

function zdot = f(z, u, parameters)
    % Flow map for Estimation over a network
A = parameters.A;
B = parameters.B;

zdot = A*z + B*u;

end

flow set C block

function inC = C(x, u, parameters) %#ok
    % Flow set indicator function.
    inC = 1;
end

For more info about Continuous-time Plant blocks, such as Continuous Process , see here .

Network:

flow map f block

function xdot = f(x, u)
    % Network flow map
    n = length(u); % measured input size
msdot = zeros(n,1); % measured continuous dynamics
jdot = 0;
tau_sdot = -1; % Timer tau_s
xdot = [msdot; jdot; tau_sdot];

end

flow set C block

function inC = C(x, vs, T_max) 
    % Flow set indicator function
    tau_s = x(end); % timer state
if tau_s >= 0 && tau_s <= T_max
    inC = 1; % report flow
elseif tau_s> T_max
    inC = 0; % do not report flow
else
    inC = 0;
end

end

jump map g block

function xplus = g(x, vs, measurement_times)
    % Jump map for network.
n = length(vs); % input size
j = x(n+1);
msplus = vs; % output = measured input
% The value tau_s is updated as a function of vs.
tau_splus = measurement_times(j+1); % Timer tau_s
j_plus = j+1;
xplus = [msplus; j_plus; tau_splus];

end

jump set D block

function inD  = D(x, vs) 
    % Jump set indicator function for network
    tau_s = x(end); % timer state
if tau_s>=0
    inD = 0; % do not report jump
elseif tau_s<= 0
    inD = 1; % report jump
else
    inD = 0;
end

end

Estimator:

flow map f block

function xdot = f(x, v, parameters)
    % Estimator flow map for Estimation over a network
    A = parameters.A;
    B = parameters.B;
    n = length(A);
    u = v(end);
    z = x(1:n);
    jdot = 0;
    zdot = A*z + B*u;
    xdot = [zdot; jdot];
end

flow set C block

function inC  = C(x, v) 
    % Estimator flow set for Estimation over a network
    j = x(end); % internal communication memory event
    jnet = v(end-2); % communication event
if j == jnet
    inC = 1; % report flow
else
    inC = 0;
end

end

jump map g block

function xplus = g(x, v, parameters)
    % Estimator jump map for Estimation over a network
    A = parameters.A;
    M = parameters.M;
    L = parameters.L;
n = length(A);
z = x(1:n);
jnet = v(end-2); % communication event
y = v(1:end-3);

jplus = jnet;
zplus = z + L*(y-M*z);
xplus = [zplus; jplus];

end

jump set D block

function inD  = D(x, v) 
    % Estimator jump set for Estimation over a network
    j = x(end); % internal communication memory event
    jnet = v(end-2); % communication event
if j ~= jnet
    inD = 1; % report jump
else
    inD = 0;
end

end

Example Output

The solution to the estimation over network system from $z(0,0)=[1, 0.1, 1, 0.6]^\top, \hat{z}(0,0)=[-10, 0.5, 0, 0, 0]^\top$ , u(t)=50*sin(0.1*t),and with T=30 , J=100 , rule=1 shows that the estimated state approaches the system's state before $t = T$ , with a small error.

References

[1] F. Ferrante, F. Gouaisbaut, R. G. Sanfelice, and S. Tarbouriech. State estimation of linear systems in the presence of sporadic measurements. Automatica, 73:101–109, November 2016.