GAMS约束编写：t≠n索引限制实现与rPower计算异常排查

核心问题排查 #

1. t≠n索引逻辑实现 #

直接在约束定义的美元控制条件中增加索引排除规则即可，不需要单独定义两个完全相同的类集合n、t，复用单个类集合就能完成索引控制，避免冗余集合带来的映射错误。

2. rPower计算结果异常原因 #

结果偏差来自三个代码错误：

• 约束未加样本-类从属控制条件：原代码对所有k、n、t组合生成约束，实际上仅当样本k属于类n时，左端的类内距离项才有意义，多余约束直接导致求解逻辑错误
• rPower数值特性问题：GAMS的rPower(x,exp)为实数幂运算，要求底数x严格非负，求解器迭代过程中一旦因为数值误差出现极小的负值，rPower会返回无意义的极大值，直接带偏迭代方向
• 幂运算未处理底数符号：原代码直接对a(k,j)-y(n,j)求p次幂，当迭代中差值为负、p为非整数时会直接返回非法值，即便p=2也可能因为数值误差触发异常

修复方案 #

p=2的场景下不需要调用rPower开根号，利用范数非负性将不等式两边做平方等价变换，完全规避幂运算的数值不稳定问题，求解效率和稳定性提升明显。
修复后可直接运行的完整代码如下：

Sets
k   sample number /1,2,3,4,5,6,7,8,9,10,11,12/
j   criterion /1,2,3,4/
n   class  /1, 2, 3/
nk(n,k) / 1.(1*4), 2.(5*8), 3.(9*12) /;

Table a(k,j)
    1   2   3   4
1   1   1   1   3
2   1   1   2   2
3   1   2   1   3
4   1   2   2   1
5   1   1   1   3
6   1   1   2   2
7   1   2   1   3
8   1   2   2   1
9   1   1   1   3
10  1   1   2   3
11  1   2   1   3
12  1   2   2   1;

Parameters 
    p /2/
    r /1e-6/
    y(n,j)
;
* 计算类中心，复用n集合无需额外定义重复的t集合
y(n,j) = sum(nk(n,k), a(k,j))/card(nk(n,k));

Variables
    e(k)
    w(j)
    z
;

Positive Variable e, w;

Equations
    constraint1
    constraint2
    constraint3
    objective
;

objective.. z =e= sum(k, e(k));
* 约束控制条件：仅对属于n类的样本k生效，排除t=n的情况
* 采用平方等价形式替代开方运算，彻底规避rPower数值异常
constraint1(n,k,t)$(nk(n,k) and ord(n)<>ord(t)).. 
    sum(j, w(j)*sqr(a(k,j)-y(n,j))) 
    =l= 
    sum(j, w(j)*sqr(a(k,j)-y(t,j))) + 2*sqrt(sum(j, w(j)*sqr(a(k,j)-y(t,j))))*(e(k)-r) + sqr(e(k)-r)
;

* 通用p范数版本（非p=2场景使用，需对差值取绝对值保证非负）
* constraint1(n,k,t)$(nk(n,k) and ord(n)<>ord(t)).. 
*   rPower(sum(j, w(j)*rPower(abs(a(k,j)-y(n,j)),p)), 1/p) - e(k)
*   =l=
*   rPower(sum(j, w(j)*rPower(abs(a(k,j)-y(t,j)),p)), 1/p) - r;

constraint2.. sum(j, w(j)) =e= 1;
constraint3(j).. w(j) =g= 1e-8;

Model nonrobust / all /;
* 换用IPOPT求解器处理光滑范数约束稳定性优于MINOS
Option nlp=IPOPT, limrow=100;
solve nonrobust minimizing z using nlp;

display z.l, e.l, w.l, y;

关键修改点说明 #

• 删除冗余的t集合、tk映射，所有类索引复用n集合，减少索引出错概率
• 给约束增加nk(n,k)控制条件，仅对类内样本生成对应约束，符合模型逻辑
• p=2场景下采用平方等价形式，完全移除rPower调用，从根源解决数值异常问题
• 通用p范数版本中对差值项加abs()取绝对值，保证rPower底数始终非负
• 调整权重下界为1e-8，避免权重取0导致的求解奇异性
• 替换NLP求解器为IPOPT，对这类光滑约束的求解精度和收敛性更好

内容的提问来源于stack exchange，提问作者diabolik