用yalmip+Gurobi求解优化模型，求解结果显示正确，但是结果不满足约束条件，不知道哪里出问题了。求帮我看看！！
代码如下：
clear
close all
clc
warning('off','YALMIP:SuspectNonSymmetry');
m=9; % 9个可行状态
E=[1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1];
O=[1 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0
0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 1];
J1=[0 0 0 1 0 0 1 0 0
0 0 0 0 1 0 0 1 0
0 0 0 0 0 1 0 0 1
0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0]; % DM1的可达矩阵
J2=[0 1 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 0 0 0
0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 1
0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0]; % DM2的可达矩阵
F11=[0.5 0.83 0.99 0.87 0.22 0.77 0.78 0.82 0.25
0.17 0.5 0.46 0.87 0.6 0.09 0.37 0.32 0.08
0.01 0.54 0.5 0.35 0.63 0.99 0.66 0.18 0.39
0.13 0.13 0.65 0.5 0.52 0.91 0.47 0.3 0.79
0.78 0.4 0.37 0.48 0.5 0.77 1 0.28 0.91
0.23 0.91 0.01 0.09 0.23 0.5 0.06 0.62 0.5
0.22 0.63 0.34 0.53 0 0.94 0.5 0.52 0.29
0.18 0.68 0.82 0.7 0.72 0.38 0.48 0.5 0.44
0.75 0.92 0.61 0.21 0.09 0.5 0.71 0.56 0.5]; % DM1中第1个相关者的原始决策矩阵
F23=[0.5 0.88 0.46 0.82 0.48 0.88 0.95 0.42 0.77
0.12 0.5 0.94 0.88 0.37 0.7 0.26 0.28 0.11
0.54 0.06 0.5 0.62 0.49 0.47 0.72 0.71 0.7
0.18 0.12 0.38 0.5 0.74 0.54 0.54 0.63 0.35
0.52 0.63 0.51 0.26 0.5 0.02 0.03 0.36 0.71
0.12 0.3 0.53 0.46 0.98 0.5 0.72 0.49 0.69
0.05 0.74 0.28 0.46 0.97 0.28 0.5 0.35 0.35
0.58 0.72 0.29 0.37 0.64 0.51 0.65 0.5 0.98
0.23 0.89 0.3 0.65 0.29 0.31 0.65 0.02 0.5]; % DM2中第3个相关者的原始决策矩阵
adjusted_F11=sdpvar(m,m,'full'); % DM1中第1个相关者的调整决策矩阵
p1=sdpvar(m,1,'full'); % DM1的调整群体意见向量
adjusted_F23=sdpvar(m,m,'full');
p2=sdpvar(m,1,'full'); % DM2的调整群体意见向量
x11=binvar(m,m,'full'); % 0-1变量，如果F11的（i，j）位置元素调整，则x11（i，j）=1，否则x11（i，j）=0
x23=binvar(m,m,'full');
P1Z=binvar(m,m,'full'); % DM1的偏好矩阵 P1+
J1J=binvar(m,m,'full'); % J1J 为 DM1 的可达改良矩阵 J1+
PF1Z=sdpvar(m,m,'full'); % DM1的偏好矩阵 P1-=
PF1ZT=sdpvar(m,m,'full'); % DM1的 P1-= 的转置
JP1=sdpvar(m,m,'full'); % JP1 为 J2*(P1-=)T
A1=binvar(m,m,'full'); % A 为 sign（J2*(P1-=)T）
B1=sdpvar(m,m,'full'); % B 为 E-sign（J2*(P1-=)T）
GMRR1=sdpvar(m,m,'full'); % DM1 的 GMR矩阵
k=128; %二进制扩展线性化相乘变量
Z1=sdpvar(m,m,k,'full');
V1=sdpvar(m,m,k,'full');
O1=sdpvar(m,m,k,'full');
Q1=sdpvar(m,m,k,'full');
y1=binvar(m,m,k,'full');
d_J1J=1/(2^k);
P2Z=binvar(m,m,'full'); % DM2的偏好矩阵 P2+
J2J=binvar(m,m,'full'); % J2J 为 DM2 的可达改良矩阵 J2+
PF2Z=sdpvar(m,m,'full'); % DM2的偏好矩阵 P2-=
PF2ZT=sdpvar(m,m,'full');% DM2的 P2-= 的转置
JP2=sdpvar(m,m,'full');% JP2 为 J1+*(P2-=)T
A2=binvar(m,m,'full');% A2 为 sign（J1+*(P2-=)T）
B2=sdpvar(m,m,'full'); % B2 为 E-sign（J1+*(P2-=)T）
SEQ2=sdpvar(m,m,'full');% DM2 的 SEQ矩阵
U2=sdpvar(m,m,k,'full');
V2=sdpvar(m,m,k,'full');
O2=sdpvar(m,m,k,'full');
M2=sdpvar(m,m,k,'full');
y2=binvar(m,m,k,'full');
d_J2J=1/(2^k);
S11=sdpvar(m,m,'full'); % DM1中第1个相关者的总调整量
S23=sdpvar(m,m,'full');
CI11=sdpvar(1,1); % adjusted_F11的一致性指数
CI23=sdpvar(1,1);
Constraints=[];
% 求解调整的群体决策矩阵和群体偏好向量
Constraints = [Constraints, sum(adjusted_F11,2)==p1]; %求解 DM1 的调整群体偏好向量
Constraints = [Constraints, sum(adjusted_F23,2)==p2]; %求解 DM2 的调整群体偏好向量
Constraints = [Constraints, p1>=0];
Constraints = [Constraints, p2>=0];
% 加性偏好关系基本约束
for i=1:9
for j=1:9
Constraints = [Constraints, adjusted_F11(i,j)+adjusted_F11(j,i)==1];
end
end
for i=1:9
for j=1:9
Constraints = [Constraints, adjusted_F23(i,j)+adjusted_F23(j,i)==1];
end
end
Constraints = [Constraints, 0<=adjusted_F11<=1];
Constraints = [Constraints, 0<=adjusted_F23<=1];
% 求解 DM1 的可达改良矩阵 J1J（J1+） DM2 的可达改良矩阵 J2J（J2+）
for i=1:9
for j=1:9
Constraints=[Constraints, -100000*(1-J1J(i,j))+0.00000001<=J1(i,j)*(p1(j)-p1(i))<=1000000*J1J(i,j)];
end
end
for i=1:9
for j=1:9
Constraints=[Constraints, -100000*(1-J2J(i,j))+0.00000001<=J2(i,j)*(p2(j)-p2(i))<=1000000*J2J(i,j)];
end
end
% 求解 DM1 的 偏好矩阵P1Z（P1+） DM2 的 偏好矩阵P2Z（P2+）
% 求解结果不满足约束条件：p1(j)-p1(i)>0时，P1Z(i,j)=1；p1(j)-p1(i)<0或=0时，P1Z(i,j)=0.该数学转化公式已验证正确
for i=1:9
for j=1:9
Constraints=[Constraints, -100000*(1-P1Z(i,j))+0.00000001<=p1(j)-p1(i)<=1000000*P1Z(i,j)];
end
end
for i=1:9
for j=1:9
Constraints=[Constraints, -100000*(1-P2Z(i,j))+0.00000001<=p2(j)-p2(i)<=1000000*P2Z(i,j)];
end
end
% 求解 DM1 的 偏好矩阵 PF1Z(P1-=) DM2 的 偏好矩阵 PF2Z(P2-=)
Constraints=[Constraints, PF1Z ==E-O-P1Z];
Constraints=[Constraints, PF2Z ==E-O-P2Z];
% 求解 DM1 的 GMR矩阵
Constraints=[Constraints, PF1ZT==transpose(PF1Z)];
Constraints=[Constraints, JP1==J2*PF1ZT];
Constraints=[Constraints, 10000*A1>=JP1];
Constraints=[Constraints, A1<=10000*JP1];
Constraints=[Constraints, B1==E-A1];
%二进制扩展求解GMRR1==J1J*B1 将J1J*B1线性化
for i=1:9
for j=1:9
for s=1:128
Constraints=[Constraints, V1(i,j,s)==B1(i,j)-Z1(i,j,s)];
end
end
end
Constraints=[Constraints, (0<=V1)&(V1<=1000*(1-y1))];
Constraints=[Constraints, (0<=Z1)&(Z1<=1000*y1)];
for i=1:9
for j=1:9
for s=1:128
Constraints=[Constraints, O1(i,j,s)==(2^(s-1))*y1(i,j,s)];
end
end
end
for i=1:9
for j=1:9
for s=1:128
Constraints=[Constraints, Q1(i,j,s)==(2^(s-1))*Z1(i,j,s)];
end
end
end
Constraints=[Constraints, J1J==(1/(2^k))*sum(O1,3)];
Constraints=[Constraints, (1/(2^k))*sum(Q1,3)==GMRR1];
% 求解 DM2 的 SEQ矩阵
Constraints=[Constraints, PF2ZT==transpose(PF2Z)];
Constraints=[Constraints, JP2==J1*PF2ZT]; % GMR 与 SEQ 差别就在这里
Constraints=[Constraints, 10000*A2>=JP2];
Constraints=[Constraints, A2<=10000*JP2];
Constraints=[Constraints, B2==E-A2];
%二进制扩展求解SEQ2==J2J*B2 将J2J*B2线性化
for i=1:9
for j=1:9
for s=1:128
Constraints=[Constraints, V2(i,j,s)==B2(i,j)-U2(i,j,s)];
end
end
end
Constraints=[Constraints, (0<=V2)&(V2<=1000*(1-y2))];
Constraints=[Constraints, (0<=U2)&(U2<=1000*y2)];
for i=1:9
for j=1:9
for s=1:128
Constraints=[Constraints, O2(i,j,s)==(2^(s-1))*y2(i,j,s)];
end
end
end
for i=1:9
for j=1:9
for s=1:128
Constraints=[Constraints, M2(i,j,s)==(2^(s-1))*U2(i,j,s)];
end
end
end
Constraints=[Constraints, J2J==(1/(2^k))*sum(O2,3)];
Constraints=[Constraints, (1/(2^k))*sum(M2,3)==SEQ2];
% 指定均衡解
Constraints=[Constraints, GMRR1(2,2)==0]; %要求第3个方案为DM1的GMRR解
Constraints=[Constraints, SEQ2(2,2)==0]; %要求第3个方案为DM2的SEQ解
% 求解每个相关者的调整量矩阵
Constraints=[Constraints, abs(F11-adjusted_F11)==S11]; % DM1中第1个相关者的调整量矩阵
Constraints=[Constraints, abs(F23-adjusted_F23)==S23];
% 求解每个相关者调整元素的个数
Constraints=[Constraints, abs(adjusted_F11-F11)<=100*x11]; %求解x11
Constraints=[Constraints, -100*(1-x11)+0.0001<=abs(adjusted_F11-F11)];%求解x11
Constraints=[Constraints, abs(adjusted_F23-F23)<=100*x23]; %求解x23
Constraints=[Constraints, -100*(1-x23)+0.0001<=abs(adjusted_F23-F23)];%求解x23
% 目标函数
z1=sum(S11(:))+sum(S23(:)); % 目标函数为最小化总调整量
% 求解优化模型
ops=sdpsettings('solver', 'gurobi', 'verbose', 0,'gurobi.NumericFocus' , 3); %求解
optimize(Constraints,z1,ops)
value(GMRR1),value(J1J),value(P1Z),value(A1),value(JP1),value(SEQ2),value(J2J),value(P2Z),value(A2),value(JP2),value(p1),value(p2)
