% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The program "Min_energy_perturbation.m" reads a square matrix, the probability transition matrix of the network, "P_matrix.xls" in excel format and returns the outputs "C", "nu", "gamma", and "mu" defined as follows: C: the minimal energy perturbation matrix, which forces the network to transition from its initial steady-state to the desired steady-state pid. nu: the Lagrange multiplier associated with the equality: pid'*(P+C) == pid' gamma: the Lagrange multiplier associated with the equality C*OneVec == ZeroVec mu: the Lagrange multiplier associated with the inequality P+C >= 0. This code uses the CVX software. For installation instructions, please refer to: http://cvxr.com/cvx/download % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0001 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0002 % The program "Min_energy_perturbation.m" reads a square matrix, the 0003 % probability transition matrix of the network, "P_matrix.xls" in excel format and returns the outputs "C", "nu", 0004 % "gamma", and "mu" defined as follows: 0005 % C: the minimal energy perturbation matrix, which forces the network to 0006 % transition from its initial steady-state to the desired steady-state pid. 0007 % nu: the Lagrange multiplier associated with the equality: pid'*(P+C) == 0008 % pid' 0009 % gamma: the Lagrange multiplier associated with the equality C*OneVec == 0010 % ZeroVec 0011 % mu: the Lagrange multiplier associated with the inequality P+C >= 0. 0012 % 0013 % This code uses the CVX software. For installation instructions, please 0014 % refer to: http://cvxr.com/cvx/download 0015 %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0016 0017 0018 %% Read the probability transition matrix "P_matrix.xls" 0019 [P,Txt,Raw]=xlsread('P_matrix'); 0020 0021 %% Global variable definition 0022 N = size(P,1); 0023 OneVec = ones(N,1); 0024 ZeroVec = zeros(N,1); 0025 0026 %% Define the desired steady-state \pi_d 0027 pid1 = 0.015525 * ones(N/2,1); 0028 pid2 = 10^(-4) * ones(N/2,1); 0029 pid = [pid1;pid2]; 0030 0031 %% The optimization algorithm using the CVX software 0032 cvx_begin 0033 variable C(N,N); 0034 dual variable nu; 0035 dual variable gamma; 0036 dual variable mu; 0037 0038 minimize( norm(C) ); 0039 subject to 0040 nu: pid'*(P+C) == pid'; % equality constraint: pid is the stationary distribution of (P+C) 0041 gamma: C*OneVec == ZeroVec; % equality constraint: C is a zero-row sum matrix 0042 mu: P+C >= 0; %inequality constraint: P+C is elementwise positive 0043 cvx_end 0044 0045 0046 0047