%% ECE 602 Project: Stochastic sensor scheduling via distributed convex optimization
% Team members: Liang Xue£¬Conghao Zhou

%% Introduction 
% In paper [1], the authors proposed a stochastic scheduling strategy for estimating the states of N discrete-time
% linear time invariant (DTLTI) dynamic systems, where only one system can be observed by the sensor
% at each time instant due to practical resource constraints. The idea of their stochastic strategy is that a
% system is randomly selected for observation at each time instant according to a pre-assigned probability
% distribution. They aim to find the optimal pre-assigned probability in order to minimize the maximal
% estimate error covariance among dynamic systems. They show that under mild conditions, the
% stochastic scheduling problem gives an upper bound on the performance of the optimal sensor selection
% problem, notoriously difficult to solve.They next relax the stochastic scheduling problem into a tractable
% suboptimal quasi-convex form. Then, they show that the new problem can be decomposed into coupled
% small convex optimization problems, and it can be solved in a distributed fashion. Finally, for scheduling
% implementation, they propose centralized and distributed deterministic scheduling strategies based on the
% optimal stochastic solution and provide simulation examples.

%% System Model and Problem formulation
% Consider that a set of $N$ DTLTI systems (targets) evolving according to the equations
%
% $x_i[k+1] = A_ix_i[k] + w_i[k]$,  $\quad\quad$  $i=1,2,...N$
%
% where $x_i[k] \in R^{n_i}$ is the process state vector and $w_i[k] \in R^{n_i}$ is assumed to be an independent Gaussian noise with zero mean
% and covariance matrix $Q_i \succ 0$. The initial state $x_i[0]$ is assumed to
% be an independent Gaussian random variable with zero mean and
% covariance matrix $\pi_i[0]$. 
%
% As a result of the sensor's limited range of sensing or the congestion of the sensing channel, at time instant $k$, only
% one system can be observed as 
%
% $$ \widetilde{y_i} = \xi_i[k](C_ix_i[k]+v_i[k])$$
%
% where $\xi_i[k]$ the indicator function indicating whether or not the
% system _i_ is observed at time instant _k_, and accordingly we have
% constraint $\Sigma_{i=1}^N \xi_i[k]=1$.  The measurement noise is $v_i[k] \in R^{p_i}$,
% which is assumed to be independent Gaussian with zero mean and
% covariance matrix $R_i \succ 0$.

%%
% Problem formulation*: Assume that  all $i \in \{1, 2, \dots,N\}$, the pair $(A_i,
% Q_i^{1/2})$ is controllable and the pair $(A_i, C_i)$ is detectable.
% Denote $\hat{x}_i[k]$ as the estimate at time $k$, obtained by a causal estimator for system $i$, which depends on the past and current observations $\{\tilde{y}_i[j]\}_{j=1}^{k}$. They began by considering problem of minimizing (in
% the limit) the maximal estimate error. The problem can be formulated mathematically as
%
% $$\min \limits_{\hat{x_i},\{\xi_i[j]\}_{j=1}^{\infty}} \max \limits_i
% \Big( \lim \sup \limits_{T \rightarrow \infty}\frac{1}{T} \sum_{k=1}^{T}
% {E}[(x_i[k]-\hat{x}_i[k])'(x_i[k]-\hat{x}_i[k])]\Big)$$
%
% $$
% \begin{array}{cl}
% s.t.  &x_i[k+1] = A_ix_i[k] + w_i[k],  \quad\quad  i=1,2,...N\\
%  &\widetilde{y_i} = \xi_i[k](C_ix_i[k]+v_i[k]),\quad\quad  i=1,2,...N\\
% &\sum_{i=1}^N \xi_i[k]=1\end{array}
% $$
%



%% Proposed Solution
% The authors relaxed the primal problem and transferred it to be a
% solvable convex optimization problem, which has the form as follows:
%
% $$\mu(\gamma)=\sum \limits_{i=1}^N \min \limits_{q_i} q_i \label{equation 1}$$
%
% subject to: $Tr(X_i(q_i)) \leq \gamma$, $q_i^c < q_i \leq 1$
%
% It is easy to infer that $\mu(\gamma)$ is deceasing w.r.t performance
% $\gamma$. Thus, the problem can be solved by using bisection on $\gamma$
% until $\mu(\gamma)=1$
% 
% The problem has _N_ independent subproblems:
%
% $$\min \limits_{q,H,G \succ 0, Y \succ 0} q$$
%
% $$
% \begin{array}{cl}
% subject to & Tr(Y) \leq \gamma, i=1, 2, \dots, N\\
% &\pmatrix{Y & I \cr I & G} \succeq \mathbf{0}\\
% &\Gamma_q(G,H) \succ 0, q^c< q \leq 1 \end{array}
% $$
%
% where $\Gamma_q(G,H)$ is
%
% $$\pmatrix{G & \sqrt{1-q}GA & \sqrt{q}(GA+HC) &\sqrt{q}H &G \cr \sqrt{1-q}GA & G &0 &0 &0 \cr \sqrt{q}(GA+HC) &0 &G &0 &0 \cr \sqrt{q}H &0 &0 &R^{-1} &0 \cr G &0 &0 &0 & Q^{-1}}$$
%
% The problem  is equivalent to a quasi-convex optimization problem that
% can be solved by using bisection on $\gamma$ until $\mu(\gamma)$ within the desired accuracy.
%
% For each level $\gamma$, each of the _N_ independent subproblems  can be solved by solving the subproblems also using bisection for variable $q_i$ and iterating _LMI_ feasibility problems.

%% Associated algorithms
% The Algorithm1 is a  distributed computing algorithm to solve the problem. The inputs of the algorithm
% are global information assumed to be known by each estimator in
% prior. Denote $\gamma^{opt}$ as the optimal objective value of the problem. The authors present the algorithm
% under the assumption that the interval [l, u] is selected to contain
% $\gamma^{opt}$. i.e., $l \leq \gamma^{opt} \leq u$ at each step. Then the algorithm is
% guaranteed to converge to the optimal objective value $\gamma^{opt}$ within
% the desired tolerance.
% 
% <<Algorithm1.PNG>>
% 
% The Algorithm2 is used to generate deterministic sequences of minimal consecutiveness. These
% sequences are periodic and switch among targets most often compatibly with the optimal scheduling distribution.
% 
% <<Algorithm2.PNG>>
% 
% The complexity of running Algorithm 2
% is O(L). In the stochastic scheduling strategy, the integer
% value $n_i$ is obtained from the optimal probability distribution,
% i.e., $n_i = \lfloor q_iL \rfloor$. In order to generate a scheduling sequence precisely
% matching the distribution $q_i's$ is within certain precision, _L_ should
% be chosen such that $q_iL(i = 1, 2,\dots N)$ is an integer. The rare
% cases where $L$ is supposed to be infinity, are approximated by taking $L$ large, and a viable sequence is then obtained by periodic
% continuation. Therefore, besides practical requirements on the
% length of measuring period, we should choose L by taking into
% account both the computing capability of the centralized scheduler
% and the probability matching precision.

%% Data Source
% In the simulation part, the authors give two examples. The first example considers a single sensor for measuring two DTLTI systems with
% the given state space representations. In example 1, the empirical performance curve
% is obtained by averaging 5000 Monte Carlo simulations. The other considers  three random-walk vehicles in an
% area and a single sensor equipped with a camera is used for tracking their 1-D positions. The dynamics of their positions are evolving as with $a_i = 1$. But they are subject to different process
% noises, measurement noises and delays. Here the authors assume that these
% vehicles have $1, 2, 2$ time-step measurement delays, respectively. 

%% Visualization of Results
% * example 1
%
% A single sensor for measuring two DTLTI systems with
% the following state space representations,
%
% $A_1=\pmatrix{0 & 1 \cr -0.49 & 1.4}$,
% $C_1=\pmatrix{1 & 0}$
%
% $Q_1=\pmatrix{5 & 0 \cr 0 & 5}$,
% $R_1=0.5$
%
% $A_2=\pmatrix{0 & 1 \cr -0.72 & 1.7}$,
% $C_1=\pmatrix{1 & 0}$
%
% $Q_1=\pmatrix{1 & 0 \cr 0 & 1}$,
% $R_1=1$ 
%
% The code to compute the $q$ is shown below.

%% 
% 
%  close all;
%  clear all;
%  A1 = [0 1;
%      -0.49 1.4];
%  C1 = [1; 0]';
%  Q1 = [5 0;
%       0 5];
%  R1 = 0.5;
%  A2 = [0 1;
%      -0.72 1.7];
%  C2 = [1; 0]';
%  Q2 = [1 0;
%       0 1];
%  R2 = 1;
% 
%  N =10;
%  epsilon = 1e-2;
%  e1 =  1e-4;
%  e2 =  1e-4;
%  q_c = 0;
%  u = 100;
%  l = 0;
%  q_u = 1;
%  q_l = q_c;
% 
%  ui = u;
%  li = l;
%  for i = 1 : 10
%     qui(1) = q_u;
%     qui(2) = q_u;
%     qli(1) = q_l;
%     qli(2) = q_l;
%     if ui - li > epsilon
%         garmma = (ui + li)/2;
%         while qui(1) - qli(1) > e1
%             q(1) = (qui(1) + qli(1))/2;
%             cvx_begin sdp quiet
%                 variable H(2,1) 
%                 variable G(2,2) 
%                 variable Y(2,2)
%                 minimize (trace(Y))
%                 subject to
%                    [Y eye(2);
%                     eye(2) G] >= 0
%                    [G, sqrt(1-q(1)).*G*A1, sqrt(q(1)).*(G*A1 + H*C1), sqrt(q(1)).*H, G;
%                            (sqrt(1-q(1)).*G*A1)', G, zeros(2,2), zeros(2,1), zeros(2,2);
%                             (sqrt(q(1)).*(G*A1+ H*C1))', zeros(2,2), G, zeros(2,1), zeros(2,2);
%                             (sqrt(q(1)).*H)', zeros(1,2), zeros(1,2), inv(R1), zeros(1,2);
%                              G', zeros(2,2), zeros(2,2), zeros(2,1), inv(Q1)] > 0
%                    G > 0
%                    Y > 0 
%             cvx_end
%             if cvx_optval <= garmma
%                 qui(1) = q(1);
%             else
%                 qli(1) = q(1);
%             end
%         end
%         q_opt(1) = q(1);
%       
%         while qui(2) - qli(2) > e2
%             q(2) = (qui(2) + qli(2))/2;
%             cvx_begin sdp quiet
%                 variable H(2,1) 
%                 variable G(2,2) 
%                 variable Y(2,2)
%                 minimize (trace(Y))
%                 subject to
%                    [G, sqrt(1-q(2)).*G*A2, sqrt(q(2)).*(G*A2 + H*C2), sqrt(q(2)).*H, G;
%                            (sqrt(1-q(2)).*G*A2)', G, zeros(2,2), zeros(2,1), zeros(2,2);
%                             (sqrt(q(2)).*(G*A2+ H*C2))', zeros(2,2), G, zeros(2,1), zeros(2,2);
%                             (sqrt(q(2)).*H)', zeros(1,2), zeros(1,2), inv(R2), zeros(1,2);
%                              G', zeros(2,2), zeros(2,2), zeros(2,1), inv(Q2)] > 0
%                    [Y eye(2);
%                     eye(2) G] >= 0
%                    G > 0
%                    Y > 0 
%             cvx_end
%             if cvx_optval <= garmma
%                 qui(2) = q(2);
%             else
%                 qli(2) = q(2);
%             end
%         end
%         q_opt(2) = q(2);
%         
%         mu = q_opt(1) + q_opt(2);
%         if mu <= 1
%             ui = garmma;
%         else 
%             li = garmma;
%         end
%     y(i) = garmma;    
%     end
%  
% end
%
%%
% By the code above, we get that $q_1 = 0.6738, q_2 = 0.3260$, which is the
% same as the result given by the authors. The optimal cost we get is
% 59.0820, which is similar to the results given in the paper 59.1. This simulation result
% verifies that result of the problem provides an effective upper bound on the expected
% estimate error covariances of the stochastic scheduling strategies.
% This is further exemplified in figure 1, which shows the empirical
% expected estimate error covariance of system.

%%
close all;
clear all;
warning off
A1 = [0 1;
     -0.49 1.4];
C1 = [1; 0]';
Q1 = [5 0;
      0 5];
R1 = 0.5;
A2 = [0 1;
     -0.72 1.7];
C2 = [1; 0]';
Q2 = [1 0;
      0 1];
R2 = 1;

epsilon = 1e-2;
e1 =  1e-4;
e2 =  1e-4;
q_c = 0;
u = 100;
l = 0;
q_u = 1;
q_l = q_c;

ui = u;
li = l;
% Calculate the optimal q_i 
for i = 1 : 10
    qui(1) = q_u;
    qui(2) = q_u;
    qli(1) = q_l;
    qli(2) = q_l;
    if ui - li > epsilon
        garmma = (ui + li)/2;
        while qui(1) - qli(1) > e1
            q(1) = (qui(1) + qli(1))/2;
            cvx_begin sdp quiet
                variable H(2,1) 
                variable G(2,2) 
                variable Y(2,2)
                minimize (trace(Y))
                subject to
                   [Y eye(2);
                    eye(2) G] >= 0
                   [G, sqrt(1-q(1)).*G*A1, sqrt(q(1)).*(G*A1 + H*C1), sqrt(q(1)).*H, G;
                           (sqrt(1-q(1)).*G*A1)', G, zeros(2,2), zeros(2,1), zeros(2,2);
                            (sqrt(q(1)).*(G*A1+ H*C1))', zeros(2,2), G, zeros(2,1), zeros(2,2);
                            (sqrt(q(1)).*H)', zeros(1,2), zeros(1,2), inv(R1), zeros(1,2);
                             G', zeros(2,2), zeros(2,2), zeros(2,1), inv(Q1)] > 0
                   G > 0
                   Y > 0 
            cvx_end
            if cvx_optval <= garmma
                qui(1) = q(1);
            else
                qli(1) = q(1);
            end
        end
        q_opt(1) = q(1);
      
        while qui(2) - qli(2) > e2
            q(2) = (qui(2) + qli(2))/2;
            cvx_begin sdp quiet
                variable H(2,1) 
                variable G(2,2) 
                variable Y(2,2)
                minimize (trace(Y))
                subject to
                   [G, sqrt(1-q(2)).*G*A2, sqrt(q(2)).*(G*A2 + H*C2), sqrt(q(2)).*H, G;
                           (sqrt(1-q(2)).*G*A2)', G, zeros(2,2), zeros(2,1), zeros(2,2);
                            (sqrt(q(2)).*(G*A2+ H*C2))', zeros(2,2), G, zeros(2,1), zeros(2,2);
                            (sqrt(q(2)).*H)', zeros(1,2), zeros(1,2), inv(R2), zeros(1,2);
                             G', zeros(2,2), zeros(2,2), zeros(2,1), inv(Q2)] > 0
                   [Y eye(2);
                    eye(2) G] >= 0
                   G > 0
                   Y > 0 
            cvx_end
            if cvx_optval <= garmma
                qui(2) = q(2);
            else
                qli(2) = q(2);
            end
        end
        q_opt(2) = q(2);
        
        mu = q_opt(1) + q_opt(2);
        if mu <= 1
            ui = garmma;
        else 
            li = garmma;
        end
    y(i) = garmma;    
    end
 
end

q1 = q(1);
q2 = q(2);
episode = 5000;
k = 50;

N = 2;
L = k;
n1 = floor(q1 * L);
n2 = floor(q2 * L);
m2 = ceil(n1/(n2 +1));
% Construct the the scheduling sequence by Algorithm 2
for i = 1 : k
    if mod(i,3) == 0
        kesei1(i) = 0;
        kesei2(i) = 1;
    else
        kesei1(i) = 1;
        kesei2(i) = 0;
    end
end

%Calulate the error covariances of empirical perfomance
%Performance of MC sequence and SS sequence
p1_SS{1} = [0,0;0,0];
p2_SS{1} = [0,0;0,0];
p1_MC{1} = [0,0;0,0];
p2_MC{1} = [0,0;0,0];
tr1_SS(1:episode,1) = trace(p1_SS{1});
tr2_SS(1:episode,1) = trace(p2_SS{1});
rand('seed',15)
for i = 1 : episode  
    for j = 1 : k  
        p1_SS{j+1} = A1*p1_SS{j}*A1' + Q1 - randsrc(1,1,[1,0;q1,1-q1]) * A1 * p1_SS{j} * C1'* inv(C1 * p1_SS{j} * C1' + R1) * C1 * p1_SS{j} * A1';
        p2_SS{j+1} = A2*p2_SS{j}*A2' + Q2 - randsrc(1,1,[1,0;q2,1-q2]) * A2 * p2_SS{j} * C2'* inv(C2 * p2_SS{j} * C2' + R2) * C2 * p2_SS{j} * A2';
        p1_MC{j+1} = A1*p1_MC{j}*A1' + Q1 - kesei1(j) * A1 * p1_MC{j} * C1'* inv(C1 * p1_MC{j} * C1' + R1) * C1 * p1_MC{j} * A1';
        p2_MC{j+1} = A2*p2_MC{j}*A2' + Q2 - kesei2(j) * A2 * p2_MC{j} * C2'* inv(C2 * p2_MC{j} * C2' + R2) * C2 * p2_MC{j} * A2';
        tr1_SS(i,j+1) = trace(p1_SS{j+1});
        tr2_SS(i,j+1) = trace(p2_SS{j+1});
        tr1_MC(i,j+1) = trace(p1_MC{j+1});
        tr2_MC(i,j+1) = trace(p2_MC{j+1});
    end
end   
tr1_emp = mean(tr1_SS);
tr1_opt(1:k) = 59.1;
%tr1_emp_cost = mean(mean(tr1_SS));
%tr1_MC_cost = mean(mean(tr1_MC));

figure(1)
plot(tr1_emp(2:k+1),'*-');
hold on
plot(tr1_opt,'k-','lineWidth',2);
axis([1 50 0 inf])
xlabel('time step (System 1)')
ylabel('Tr(P1(k))')
legend('Empirical performance','Upper bound')
snapnow
figure(2)
subplot(2,1,1)
plot(tr1_SS(1,:),'b','linewidth',2);
axis([1 50 0 inf])
hold on;
plot(tr1_MC(1,:),'r--','linewidth',2);
axis([1 50 0 inf])
xlabel('time k (System 1)')
ylabel('Tr(P1(k))')
legend('SS sequence','MC sequence')

subplot(2,1,2)
plot(tr2_SS(1,:),'b','linewidth',2);
axis([1 50 0 inf])
hold on;
plot(tr2_MC(1,:),'r--','linewidth',2);
axis([1 50 0 inf])
xlabel('time k (System 2)')
ylabel('Tr(P2(k))')
legend('SS sequence','MC sequence')
snapnow

%%
% Figure 1 shows the empirical estimate error covariances for DTLTI system. The empirical performance curve
% is obtained by averaging 5000 Monte Carlo simulations. The figure we
% generated is the same as the figure given by the authors. 
%
% Figure 2 shows the  performance of minimal consecutiveness sequence vs. stochastic scheduling
% sequence. In figure 2, the blue solid curves show the evolution
% of the error covariance w.r.t time k. We see several high peaks as
% a result of consecutive loss of observations. The dash lines show
% that the MC sequence provides much smoother error covariance
% which is quite necessary in certain applications. Furthermore, the
% objective value provided by MC sequence, 55.7, is smaller than the
% empirical expected estimate error covariances, 58.7, provided by
% the SS sequences. This performance enhancement is the result of
% "smoothness" provided by the MC sequence. In a word, the MC
% sequence provides a better upper bound than that provided by the
% SS sequences.
%
% * example 2
%
% Consider three random-walk vehicles in an area and a single sensor equipped with a camera is used for tracking their $1-D$ positions. The dynamics of their positions are evolving with ai = 1. But they are subject to different process
% noises, measurement noises and delays. Here these vehicles have 1, 2 and 2 time-step measurement delays, respectively. The expanded state space representations are as follows.
%
% $A_1=\pmatrix{0 & 1 \cr 0 & 1}$,
% $B_1=\pmatrix{0 \cr 1}$
% $C_1=\pmatrix{1 & 0}$
%
% $A_i=\pmatrix{0 & 1 & 0 \cr 0 & 0 &1 \cr 0 &0 &1}$,
% $B_i=\pmatrix{0 \cr 0 \cr 1}$
% $C_i=\pmatrix{1 & 0 & 0}$
% $i=2,3$
%
% $Q_1=1$, $Q_2=2$, $Q_3=5$
% $R_1=R_2=R_3=1$.

close all;
clear all;
warning off
A1 = [0 1;
      0 1];
B1 = [0; 1];
C1 = [1; 0]';
Q1 = 1;
R1 = 0.5;
A2 = [0 1 0;
      0 0 1; 
      0 0 1];
B2 = [0; 0; 1];
C2 = [1; 0; 0]';
Q2 = 2;
R2 = 1;
A3 = [0 1 0;
      0 0 1; 
      0 0 1];
B3 = [0; 0; 1];
C3 = [1; 0; 0]';
Q3 = 2;
R3 = 1;
q1 = 0.0649;
q2 = 0.1612;
q3 = 0.7739;

k = 300;
a = 1;
rand('seed',2)
X1(:,1) = [0;40];
X2(:,1) = [0; 0; 0];
X3(:,1) = [0; 0; -10];
y1(1) =  randsrc(1,1,[1,0;q1,1-q1]) * (C1 * X1(:, 1) + normrnd(0,R1));
y2(1) =  randsrc(1,1,[1,0;q2,1-q2]) * (C2 * X2(:, 1) + normrnd(0,R2));
y3(1) =  randsrc(1,1,[1,0;q3,1-q3]) * (C3 * X3(:, 1) + normrnd(0,R3));
for i = 1 : k
    X1(:, i+1) = A1 * X1(:, i) + B1 * normrnd(0,Q1);
    X2(:, i+1) = A2 * X2(:, i) + B2 * normrnd(0,Q2);
    X3(:, i+1) = A3 * X3(:, i) + B3 * normrnd(0,Q3);
    kesei1 = randsrc(1,1,[1,0;q1,1-q1]);
    kesei2 = randsrc(1,1,[1,0;q2,1-q2]);
    kesei3 = randsrc(1,1,[1,0;q3,1-q3]);
    if kesei1 == 0
        y1(i+1) = y1(i);
    else
        y1(i+1) = kesei1 * (C1 * X1(:, i+1) + normrnd(0,R1));
    end         
    if kesei2 == 0
        y2(i+1) = y2(i);
    else
        y2(i+1) = kesei2 * (C2 * X2(:, i+1) + normrnd(0,R2)); 
    end   
    if kesei3 == 0
        y3(i+1) = y3(i);
    else
        y3(i+1) = kesei3 * (C3 * X3(:, i+1) + normrnd(0,R3));
    end
end

figure(1)
plot(X1(2,6:k+1),'k','linewidth',2);
hold on;
plot(y1(6:k+1),'r','linewidth',2);
hold on;

plot(X2(2,5:k+1),'k','linewidth',2);
hold on;
plot(y2(6:k+1),'r','linewidth',2);
hold on;

plot(X3(2,6:k+1),'k','linewidth',2);
hold on;
plot(y3(6:k+1),'r','linewidth',2);
hold on;
axis([0 300 -inf inf])
xlabel('time k')
ylabel('Position')
snapnow;
%     snapnow;
%
%%
% Figure 3 shows the positions of random-walk vehicles and tracking paths. The flat segments
% imply that no measurement is taken in the time slot and the estimator simply
% propagates the state estimate of the previous time-step. Figure 3 we generated is slightly different from the one given in the article. We think the reason might be 
% the positions of the random walk vehicles are generated  randomly. So there is some randomness here. By running the proposed distributed algorithms, the solutions of
% OP is $[q1, q2, q3] = [0.0649, 0.1612, 0.7739]$. We assume the
% existence of a centralized scheduler which constructs a random
% scheduling sequence accordingly. The tracking paths (red curves) are shown in
% comparison with actual time-varying positions (black curves) of
% random-walk vehicles. The flat segments of red curves imply that
% no measurement is taken in this time slot and the estimator simply
% propagates the state estimate of the previous time-step. In figure 3,
% it is shown that the actual path of vehicle 1 changes slowly, correspondingly the attention given to this sensor is small and the estimate path is updated rarely. In comparison, the tracking paths
% of vehicle 2 and vehicle 3 match the actual positions much better
% even though the flat segments in tracking path of vehicle 2 is occasionally visible.

%% Result Analysis & Conclusion
% The results exactly validate the results in paper [1]. we reproduce the
% results reported in the paper and the algorithms behave as described by
% the authors. The conclusions we obtained support that made by authors. In
% example 1, The simulation result
% verifies that result of the problem provides an effective upper bound on the expected
% estimate error covariances of the stochastic scheduling strategies. The
% MC sequence provides much smoother error covariance, the performance enhancement is the result of
% smoothness provided by the MC sequence. In example 3, it is shown that the actual path of vehicle 1 changes slowly, correspondingly the attention given to this sensor is small and the estimate path is updated rarely. In comparison, the tracking paths
% of vehicle 2 and vehicle 3 match the actual positions much better.
%
% In conclusion, this paper proposed a stochastic sensor scheduling
% problem where the states of N targets (i.e., DTLTI systems) need
% to be estimated. Due to some practical sensing constraints, at each
% time instant, only one target can be measured by the sensor. The
% basic idea of stochastic scheduling is that at each time instant
% a target is randomly chosen for observation according to a preassigned probability distribution. The authorrs relaxed the problem
% to a convex optimization problem and proposed a distributed
% algorithm to obtain the optimal probability distribution with the
% objective to minimize the maximal estimate error among all
% targets. 

%% Reference
% [1] C. Li and N. Elia. "Stochastic sensor scheduling via distributed convex optimization." Automatica, 58:173-182, 2015.
