%% Malicious Data Attacks on Smart Grid State Estimation: Attack Strategies and Countermeasures
%
%% Problem formulation
%
% The modern smart grid depends heavily on communication, with the purpose
% of making the power system more efficient and flexible. However, this
% integration with the cyber infrastrucutre makes the network more
% vulnerable to cyber attacks. The authors of this paper investigate the
% problem of constructing malicious data attack of smart grid state
% estimation with countermeasures that detect the presence of such attacks.
% The attacks are constructed such that to minimize the residue energy of
% attack while guaranteeing a certain level of increase of mean square
% error. For attacks detection, the Generalized Likelihood Ratio Test
% (GLRT) is used, and regularized by an $L_1$ norm penalty of the strength
% of the attack. 
%% Proposed solution
% The first step is to represent the power system, achieved by adopting a graph theoretic model. 
% The power system is composed of Buses, transmission lines, and power flow meters. 
% The goal of state estimation is to recover the full state of system, i.e
% in case of power system the voltage and phase and of every bus in the network. 
% Adopting a DC power flow model, which is a linearized version of the AC
% power flow model:
% 
% <<DCpower.png>>
% 
% where $z \in \mathbf{R}^m$ is the vector power flow measurements, $\mathbf{x} \in \mathbf{R}^n$
% the system state (voltage phase angles), $\mathbf{e}$ the Gaussian measurement noise with zero mean
% and covariance matrix $\Sigma_e$, and  vector $\mathbf{a}$ is malicious data injected by adversary.
%
% The authors consider a Bayesian framework, with random state variables
% with Gaussian distribution $\mathcal{N} (\mu_x,\Sigma_x)$, with $\mu_x=0$. In the case there is no attack 
% i.e $\mathbf{a=0}$, The minimum mean square error (*MMSE*) estimator of the sate vector $\mathbf{x}$ is a linear
% estimator given by:
% 
% <<xhat.png>>
% 
% 
% where
% 
% <<k.png>>
%
%
% The minimum mean square error, in the absence of attack is given by:
% 
% <<E0.png>>
% 
% In the case of presence of an attack, then the estimator
% $\mathbf{\hat{x} = Kz}$ is a naive estimator that ignores the possibility
% of attack. The mean square error (*MSE*) in this case is given by:
% 
% <<MSE.png>>
% 
% where the impact on the estimator from a particular attack $\mathbf{a}$ is
% the term $\Vert \mathbf{Ka} \Vert_2^2$. 
%
% For the detection problem at the control center, the authors assume a
% Bayesian model of starte variables with a multivariable Gaussian
% distribution $\mathbf{e \sim  \mathcal{N} (0,\Sigma_e)}$. 
%
% The following composite binary hypothesis are considered under the
% observation model of the DC power flow:
% 
% <<H0H1.png>>
% 
% The goal is to design a detector $\delta$: $\mathbf{R} \Longrightarrow
% \lbrace 0,1 \rbrace$, with $\delta(\mathbf{z})=1$ indicating a detection of attack
% $\mathcal({H}_1)$ and $\delta(\mathbf{(z)}=0$ the null hypothesis $\mathcal({H}_0)$
%
% 
% * *Generalized Likelihood Ratio Detector (GLRT) with $L_1$ Norm
% Regularization*
%
% There exits no uniformly most powerful test for the given composite
% binary hypotheses, therefore the authors propose a generalized likelihood
% ratio test (GLRT). The distribution of the measurement $\mathbf{z}$ under
% the two hypotheses differ only in their means:
% 
% <<distribution.png>>
% 
% where  $\mathbf{\Sigma_z} \equiv \mathbf{H \Sigma_x H^T + \Sigma_e}$. From that the GLRT can be written as:
%
% <<decision.PNG>>
% 
%
% Thus the GLRT reduces to solving
% 
% <<minimize.PNG>>
% 
%
% The $L_0$ in the above optimization problem makes it non-convex, but for
% a small k it can be solved exactly by exhaustively searching through all
% sparsity patterns. For larger k, this is not feasible, however it is a
% well known that the above optimization problem can be approximated by the
% following convex optimization
% 
% <<relaxation.PNG>>
%
%
% where the $L_1$ norm constraint is a heuristic replacement for the
% sparsity of $\mathbf{a}$. By solving the above optimization problem
% multiple times, the solution will involve an $\mathbf{a}$ with a sparsity
% $k$ with adjustment of the constant $\nu$.
%% Data sources 
% 
% The DC power flow model data used by the authors is obtained by using a
% DC State Estimation Model, using the Weighted Least Squares Method (WLS).
% The DC model is a linearized version of the AC model. The assumptions in
% this model are that Buses voltages are known and equal to 1.0 pu, also
% all shunt elements are branch resistances are neglected. Power flow is
% approximated by:
% 
% $$P_{ik}=\frac{\theta_i - \theta_k}{x_{ik}}+ e$$
% 
% where $x_{ik}$ is the reactance of branch $i-k$, $\theta_i$ is the phase
% angle at bus i, and $e$ is the measurement error. Also, power injection
% measurement at bus $i$ is given by:
% 
% $$P_i=\sum_{j \in N_j} P_{ij} +e$$
%
% where $N_j$ is the set of buses connected to bus j. Hence, the DC model for the real power measurements can be expressed
% in matrix form as follows:
% 
% $$z =  H \theta + e$$
% 
% where $z$ includes flow and injection measurements, $H$ is a function of
% the branch reactances only, and $e$ is the vector of random errors. The
% code used to generate the power flow measurements is obtained by
% following the algorithm described in [1].
%% Solution
clear;clc;
load Power_DC

%GLRT detector
%1- Sparsity =2, SNR=10db
k=2; SNR=10;
[a_1,Co_z1,Co_x1,Co_e1] = GLRT(k,z,SNR,H,x);

% %2- Sparsity =3, SNR=10db
% k=3; SNR=10;
% [a_2,Co_z2,Co_x2,Co_e2] = GLRT(k,z,SNR,H,x);
% 
% %3- Sparsity =3, SNR=6db
% k=3; SNR=6;
% [a_3,Co_z3,Co_x3,Co_e3] = GLRT(k,z,SNR,H,x);
% 
% %4- Sparsity =6, SNR=10db
% k=6;SNR=10;
% [a_4,Co_z4,Co_x4,Co_e4] = GLRT(k,z,SNR,H,x);
%% Visualization of results
% Two types of plots are used to show the performance of the detector. The
% first is the Receiver Operating Characteristics (ROC) that characterize
% the tradeoff between the probability of attack detection vs. the
% probability of false alarm. In addition, the Attacker Operating
% Characteristics (AOC)is plotted. The AOC charaterizes the tradeoff
% between the probability of being detected vs. resulting (extra)
% mean-square error at the state estimator (for a fixed probability of
% false alarm)

%1- ROC & AOC of GLRT for 2 sparsity case
%ROC: MSE with attack is 8db,SNR=10db
%AOC: Flase alarm rate=0.05, SNR=10db

k=2;SNR=10; ROC(a_1,Co_z1,Co_x1,Co_e1,k,SNR,z,H); % Sparsity k=2, SNR=10
% k=3;SNR=10; ROC(a_2,Co_z1,Co_x2,Co_e2,k,SNR,z,H); % Sparsity k=3, SNR=10
% k=3;SNR=6 ; ROC(a_3,Co_z1,Co_x3,Co_e3,k,SNR,z,H); % Sparsity k=3, SNR=6
% k=6;SNR=10; ROC(a_4,Co_z1,Co_x4,Co_e4,k,SNR,z,H); % Sparsity k=6, SNR=10


%% Analysis and conclusions
%
% 
% * A dectector based on the GLRT is implemented to detect the attack
% vector a. The optimization problem has a constraint of L_0 norm, making
% it a non convex, therefore a relaxation is achieved by approximation by
% using L_1 norm. 
%
% * The authors used L_0 and used exhaustive search over all possible
% sparsity patterns, which is only feasible for low sparsity patterns and
% the relaxation of L_1 in higher sparsity. 
% However, our implementation used only the L_1 norm, for all cases. This
% solution is feasible for all test cases, but will a more approximated
% results of low sparsity (k=2, k=3) that the results of the paper. This is why 
% it was required to solve the optimization problem several times adjusting
% the constant of L_1 to achieve the required sparsity.
%
% * The plot of ROC don't match exactly the results of the paper.
% We believe the implementation of the optimization problem is correct,
% and the problem lies in plotting the ROC curves. 
%
%% A list of source files/functions necessary for algorithm implementation
% 1- DC power state estimation data is obtained from: WLS.m in Power Flow
%  folder, the results are stored in Power_DC.mat to be used
%
% 2- main.m uses two functions:
%  1- GLRT.m: contains optimization problem code
%  for the given case of sparsity and SNR
%  2- ROC.m: plot the ROC and AOC curves for the
%  given case of sparsity and SNR
%% References
%
% [1] Abur, Ali, and Antonio Gomez Exposito. Power system state estimation: theory and implementation. CRC press, 2004.
% [2] Kosut, Oliver, et al. "Malicious data attacks on the smart grid." IEEE Transactions on Smart Grid 2.4 (2011): 645-658.

