编程语言
首页 > 编程语言> > SMT算法运行记录

SMT算法运行记录

作者:互联网

仿真模型,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为:[41223777.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,

标签:counter,记录,SMT,残差,slackIndex,算法,传感器,active,constraints
来源: https://www.cnblogs.com/creationMarvel/p/15933764.html