仿真模型，p=7,s=3;n=4。

A=[0.9966 0.0100 0.0012 0.0000;
　　-0.6730 0.9966 0.2480 0.0012;
　　0.0033 0.0000 1.0028 0.0100;
　　0.6525 0.0033 0.5578 1.0028;];
B=[0.0022 0.4475 -0.0002 -0.0461]';
C=[1 0 0 0 1 0 1;
　　0 1 0 0 0 2 0;
　　0 0 1 0 2 0 0;
　　0 0 0 1 0 1 1;]';

系统输入

u = K*x+5*cos(0.002*t);

初始状态，也是我们算法的目标结果：

x0 = [5 0 -5 0]';

受到攻击的传感器为1、2、5号，攻击信号设为[-50,50]范围内服从正态分布的随机数：

attacked_sensor_index = [1,2,5];

y(attacked_sensor_index) = y(attacked_sensor_index) + 50*randn(1) ;

收集连续tau=n=4个时刻的测量数据，减去已知输入的影响后得到如下结果，每一行代表一个时刻，每一列代表一个传感器：

我们需要根据这4个时刻、7个传感器的测量数据，估计出初始时刻的状态。下面的算法1是本算法的主体架构：

对应算法1、2行的初始化步骤如下：

xhat = [];
sensorsUnderAttack = [];
smt.agreeableClauses = {};
smt.conflictClauses = {};
solved = 0;

第一次运行SAT求解器，此时的约束只有“攻击数目不超过s_”。

convexConstraintAssignment  = smt.SATsolver.model;  %  获取一个可能的分配

上述指令的运行结果为[-1,-2,-3,-4,-5,-6,-7]。这里需要说明的是：结果中正元素代表我们假设该传感器受到攻击，负元素表示我们假设该传感器安全。显然，该结果含义是7个传感器均安全，相当于算法1中第四行求出的b向量为[0,0,0,0,0,0,0]，这倒也符合此时的约束，但肯定是错误的。

下一步，我们要对“我们假设的安全传感器组”进行检验，不难看出，下面两条指令的结果分别是[1,2,3,4,5,6,7]和[]。

constraints = find(convexConstraintAssignment < 0);                      % 假设未受攻击的传感器，即I集合
sensorsUnderAttack = find(convexConstraintAssignment > 0);       % 假设受到攻击的传感器

constraints这个变量就是我们送到算法2中进行验证的集合，即算法2的输入参数I。算法2如下：

 

 

 首先计算得到YI和OI，然后根据最小二乘计算x的估计值，由于此时使用的数据包含攻击，所以结果必然不准确：

Y_active = [];
O_active = [];
for counter = 1 : length(constraints)
　　Y_active = [Y_active; Y{constraints(counter)}];
　　O_active = [O_active; O{constraints(counter)}];
end
xhat = pinv(O_active)*Y_active;       % 直接求解x    结果为 [-1.85;-14.58;-14.8;24.39]，与目标值[5 0 -5 0]相差甚远。

计算各个传感器测量残差的二范数（后面进行排序），相加后即为总残差，随后的判断逻辑与伪代码中类似（考虑了噪声）：

for counter = 1 : length(constraints)
　　slack(constraints(counter)) = (norm(Y{constraints(counter)} - O{constraints(counter)}*xhat))^2; % 这里为什么没有平方？
end
if sum(slack) <= tolerance + sum(noiseBound(constraints))
　　solved = 1; % 置1 说明问题已解决  
　　return;
end

计算结果分残差向量slack为：[4122；3777.8；41.5239；919.1420；3633.2；382.9694；523.488]。可以发现，三个受攻击的传感器对应的残差明显较大。

总残差为19523，后面的阈值为0.7，说明问题没有被解决，即算法1中的status仍然是UNSAT，于是下一步我们根据刚刚计算的残差来获取新证书，即T-CERTIFICATE算法：

 

 

 对constraints变量中的传感器对应的残差从小到大排序，由于第一次循环constraints中包含了所有传感器，所以可以直接对slack变量中的元素排序，得到的slackIndex为：3,6,7,4,5,2,1。

% 1/ 基于残差排序
[~, slackIndex] = sort(slack(constraints), 'ascend');
slackIndex = constraints(slackIndex); % slackIndex为I中传感器序号的集合 按照残差升序排列

为什么这里不能直接假定残差最大的三个传感器受到攻击？

然后找到残差最小的p-2s个（准确的说应该是p-s-s^个，正常情况下传入的I中应该只包含p-s个元素，但本次循环是第一次所以例外）传感器和残差最大的传感器，这里显然分别为[3,6,7,4]和[1]。

% 2/ 找到残差最小的p-2s个传感器
indexLowSlackSensors = slackIndex(1 : end - s);          % 找到slackIndex中的前p-2s个
[~, indexSortedLowSlackSensors] = sort(dimNull(indexLowSlackSensors), 'ascend');    % 计算一下各自的零空间，维数越低价值越高
indexSortedLowSlackSensors = indexLowSlackSensors(indexSortedLowSlackSensors);
max_slack_index = slackIndex(end);    % 找到slackIndex中的最后一个

让着p-2s+1个传感器组成一个新集合I'，用算法2做验证。由于它们残差差距非常大，所以大概率结果是UNSAT，这样我们就得到了一个“冲突性证书”。

if(norm(Y_active - Y_hat) > tolerance + sum(noiseBound(indexSortedLowSlackSensors(1:counter))) )
　　% Conflict discovered
　　conflicts = [max_slack_index, indexSortedLowSlackSensors(1:counter)'];
　　foundConflict = 1;
　　break;
end

上面的判断条件中，norm(Y_active - Y_hat)的值为50.9944，