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为:[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,
标签:counter,记录,SMT,残差,slackIndex,算法,传感器,active,constraints 来源: https://www.cnblogs.com/creationMarvel/p/15933764.html