[学习笔记]Lingo超经典案例大全

LINGO是Linear Interactive and General Optimizer的缩写，即“交互式的线性和通用优化求解器”。Lingo超强的优化计算能力在很多方面（线性规划、非线性规划、线性整数规划、非线性整数规划、非线性混合规划、二次规划等）比matlab、maple等强得多，Lingo编程简洁明了，数学模型不用做大的改动（或者不用改动）便可以直接采用Lingo语言编程，十分直观。

Lingo模型由4个段构成：

（1）集合段（sets  endsets）；（2）数据段（data  enddata)；

 (3)初始段（init endinit）；（4）目标与约束段。

Lingo的五大优点：

1. 对大规模数学规划，LINGO语言所建模型较简洁，语句不多；

2. 模型易于扩展，因为@FOR、@SUM等语句并没有指定循环或求和的上下限，如果在集合定义部分增加集合成员的个数，则循环或求和自然扩展，不需要改动目标函数和约束条件；

3. 数据初始化部分与其它部分语句分开，对同一模型用不同数据来计算时，只需改动数据部分即可，其它语句不变；

4. “集合”是LINGO有特色的概念，它把实际问题中的事物与数学变量及常量联系起来，是实际问题到数学量的抽象，它比C语言中的数组用途更为广泛。 

5. 使用了集合以及@FOR、@SUM等集合操作函数以后可以用简洁的语句表达出常见的规划模型中的目标函数和约束条件，即使模型有大量决策变量和大量数据，组成模型的语句并不随之增加．

 一、求解线性整数规划、非线性整数规划问题：

1.线性整数规划：

model:

max=x1+x2;

x1+9/14*x2<=51/14;

-2*x1+x2<=1/3;

@gin(x1);@gin(x2);

end

求得x1=3，x2=1，最大值为4.运用matlab求时可以发现有两组解：x1=3，x2=1和x1=2，x2=2。通过验证也可知这两组解均满足。Lingo的一个缺陷是：每次只能输出最优解中的一个（有时不只一个）。那么，怎样求得其他解呢？一个办法是将求得的解作为约束条件，约束x1不等于3，x2不等于1，再求解。如下：

model:

max=x1+x2;

x1+9/14*x2<=51/14;

-2*x1+x2<=1/3;

@gin(x1);@gin(x2);

@abs(x1-3)>0.001;

@abs(x2-1)>0.001;

end

求得x1=2，x2=2.若再次排除这组解，发现Lingo解不出第三组解了，这时我们可以断定：此优化模型有两组解：

x1=3，x2=1和x1=2，x2=2.

求解模型时需注意：Lingo中，默认变量均为非负；输出的解可能是最优解中的一组，要判断、检验是否还有其他解（根据具体问题的解的情况或用排除已知最优解的约束条件法）。

2、非线性整数规划：

model:

sets:

row/1..4/:b;

col/1..5/:c1,c2,x;

link(row,col):a;

endsets

data:

c1=1,1,3,4,2;

c2=-8,-2,-3,-1,-2;

a=1 1 1 1 1

1 2 2 1 6

2 1 6 0 0

0 0 1 1 5;

b=400,800,200,200;

enddata     

max=@sum(col:c1*x^2+c2*x);

@for(row(i):@sum(col(j):a(i,j)*x(j))

@for(col:@gin(x));

@for(col:@bnd(0,x,99));

End

求得：x1=50,x2=99,x3=0,x4=99,x5=20.最大值为51568。

这里,我们看不出是否还有其他解，需要将已知的最优解排除掉。

利用1的方法分别可得到其他解：

x1=48,x2=98,x3=1,x4=98,x5=19.最大值为50330。

x1=45,x2=97,x3=2,x4=97,x5=18.最大值为49037。

x1=43,x2=96,x3=3,x4=96,x5=17.最大值为47859。

x1=40,x2=95,x3=4,x4=95,x5=16.最大值为46636。

......

发现x1，x2，x4，x5均单调减少，x3单调增加。最大值越来越小。可以简单判断第一组为最优的。当然，能够一一检验最好。

 

二、最优选择问题

某钻井队要从10个可供选择的井位中确定5个钻井探油，使总的钻探费用为最小。若10个井位的代号为s1,s2,...,s10,相应的钻探费用c1,c2,...,c10为5,8,10,6,9,5,7,6,10,8.并且井位选择上要满足下列限制条件：
（1） 或选择s1和s7,或选择钻探s9；
（2） 选择了s3或s4就不能选s5，或反过来也一样；
（3） 在s5,s6,s7,s8中最多只能选两个.

试建立这个问题的整数规划模型,确定选择的井位。

取0-1变量s_i，若s_i=1,则表示选取第i个井，若s_i=0,则表示不选取第i个井。建立数学模型如下：

model:
sets:
variables/1..10/:s,cost;
endsets
data:
cost=5 8 10 6 9 5 7 6 10 8;
enddata
min=@sum(variables:cost*s);
(s(1)+s(7)-2)*(s(9)-1)=0;
s(3)*s(5)+s(4)*s(5)=0;
@sum(variables(i)|i#ge#5#and#i#le#8:s(i))<=2;
@sum(variables:s)=5;
@for(variables:@bin(s));
end

求得:

                          Total solver iterations:   26

                       Variable           Value        Reduced Cost
                          S( 1)        1.000000           -4.000000
                          S( 2)        1.000000            0.000000
                          S( 3)        0.000000            2.000000
                          S( 4)        1.000000           -2.000000
                          S( 5)        0.000000            0.000000
                          S( 6)        1.000000           -1.000000
                          S( 7)        1.000000            0.000000
                          S( 8)        0.000000            0.000000
                          S( 9)        0.000000            2.000000
                         S( 10)        0.000000            0.000000

                         Objective value:   31.00000

即选择井S1，S2，S4，S6，S7以达到最小费用31.

 

 三、路径和最短问题：

设平面上有N个点，求一点，使得这个点到所有点距离之和最小。这里，取N=8。数据点是1~5的随机数。

Lingo：

model:
sets:
position/1..8/:x,y;
ab/1/:a,b;
endsets
data:
@text('E:\matlab7.0\work\data.txt')=x,y;!读入到matlab的工作空间中;
@text('E:\matlab7.0\work\data1.txt')=a,b;
enddata
x(1)=1+4*@rand(0.12345);
y(1)=1+4*@rand(0.25);
@for(position(i)|i#ge#2:x(i)=1+4*@rand(x(i-1)));!随机产生1~5中的8个点;
@for(position(i)|i#ge#2:y(i)=1+4*@rand(y(i-1)));
[obj]min=@sum(position(i):@sqrt((x(i)-a(1))^2+(y(i)-b(1))^2));!目标函数;
@bnd(1,a(1),5);
@bnd(1,b(1),5);
end

matlab：

clear;
clc;
close all;
load('data.txt');
load('data1.txt');
hold on;
plot(data1(1),data1(2),'o','MarkerSize',15,'MarkerFaceColor','r');
plot(data(:,1),data(:,2),'or','MarkerSize',15,'MarkerFaceColor','b');
set(gcf,'Color','w');
set(gca,'FontSize',16)
grid off;
data1=repmat(data1,8,1);
P=[data1(:,1)';data(:,1)'];
Q=[data1(:,2)';data(:,2)'];
plot(P,Q,'g','LineWidth',2);
xlabel('x');
ylabel('y');
title('Solving the problem of the minimun distance of tne sum of all the blue points towards the being known red point.');
gtext(['The minimun distance is ',num2str(10.2685),'.'],'FontSize',16,'Color','r');

 

三、运输+选址问题：

 某公司有6个建筑工地，位置坐标为(ai, bi) (单位：公里),水泥日用量di (单位：吨）

i        1         2          3          4          5           6

a       1.25      8.75       0.5        5.75        3           7.25

b       1.25      0.75       4.75       5           6.5         7.75

d       3          5          4         7           6            11

(1)现有2料场，位于A (5, 1), B (2, 7),记(xj,yj),j=1,2, 日储量ej各有20吨。

假设料场和工地之间有直线道路，制定每天的供应计划，即从A, B两料场分别向各工地运送多少吨水泥，使总的吨公里数最小。

取决策变量c_ij表示i工地从j料场运来的水泥量。模型（线性模型）为：

model:
sets:
demand/1..6/:a,b,d;
supply/1..2/:x,y,e;
link(demand,supply):c;
endsets
data:
a=1.25 8.75 0.5 5.75 3 7.25;
b=1.25 0.75 4.75 5 6.5 7.75;
d=3 5 4 7 6 11;
x=5 2;
y=1 7;
e=20 20;
enddata
[obj]min=@sum(link(i,j):c(i,j)*@sqrt((a(i)-x(j))^2+(b(i)-y(j))^2));!目标函数;
@for(demand(i):@sum(supply(j):c(i,j))=d(i));
@for(supply(j):@sum(demand(i):c(i,j))<=e(j));
end

求得：

C( 1, 1)        3.000000            

C( 1, 2)        0.000000            

C( 2, 1)        5.000000           

C( 2, 2)        0.000000           

C( 3, 1)        0.000000          

C( 3, 2)        4.000000                                  

C( 4, 1)        7.000000                            

C( 4, 2)        0.000000                                

C( 5, 1)        0.000000                                 

C( 5, 2)        6.000000                                 

C( 6, 1)        1.000000                             

C( 6, 2)        10.00000           
Objective value:    136.2275

（2） 改建两个新料场，需要确定新料场位置(xj,yj)和运量cij ，在其它条件不变下使总吨公里数最小。

模型一样，未知量变为料场位置(xj,yj)和运量cij ，变为非线性优化问题。

model:
sets:
demand/1..6/:a,b,d;
supply/1..2/:x,y,e;
link(demand,supply):c;
endsets
data:
a=1.25 8.75 0.5 5.75 3 7.25;
b=1.25 0.75 4.75 5 6.5 7.75;
d=3 5 4 7 6 11;
e=20 20;
enddata
init:
x=5 2;
y=1 7;
endinit
[obj]min=@sum(link(i,j):c(i,j)*@sqrt((a(i)-x(j))^2+(b(i)-y(j))^2));!目标函数;
@for(demand(i):@sum(supply(j):c(i,j))=d(i));
@for(supply(j):@sum(demand(i):c(i,j))<=e(j));
@for(supply:@free(x);@free(y));
end

求得：

C( 1, 1)        3.000000         

C( 1, 2)        0.000000                                  

C( 2, 1)        0.000000                                 

C( 2, 2)        5.000000                                  

C( 3, 1)        4.000000                                  

C( 3, 2)        0.000000                                  

C( 4, 1)        7.000000                                 

C( 4, 2)        0.000000                                   

C( 5, 1)        6.000000                                 

C( 5, 2)        0.000000                                  

C( 6, 1)        0.000000                                  

C( 6, 2)        11.00000           

（x1，y1）=（3.254884，5.652331） 

（x2，y2）=（7.250000，7.750000）

Objective value:   85.26604

四、路径最短问题：

如上图，求从S到T的最短路径。设d(x,y):城市x与城市y之间的直线距离;L(x):城市S到城市x的最优行驶路线的路长。模型为：

min {L(x)+d(x,y)}

L(S)=0

 

model:
sets:
city/S,A1,A2,A3,B1,B2,C1,C2,T/:L;
road(city,city)/S,A1 S,A2 S,A3 A1,B1 A1,B2 A2,B1 A2,B2 A3,B1 A3,B2 B1,C1 B1,C2 B2,C1 B2,C2 C1,T C2,T/:d;
endsets
data:
d=6 3 3
6 5 8 6 7 4
6 7 8 9
5 6;
L=0,6,3,3,,,,,;
enddata
@for(city(j)|j#gt#@index(city,S):L(j)=@min(road(i,j):L(i)+d(i,j)));
end

求得最短路径为20.

 

五、指派问题(0-1规划问题）：

四个人完成4项任务所用的时间如下，问如何指派任务使得完成所有任务的时间最短？

     任务   t1    t2    t3    t4

人员

m1          2     15    13    4

 

m2          10     4    14    15    

 

m3          9      14    16    13

 

m4          7      8     11     9

 

c_ij:表示第i个人完成第j项任务所用的时间；

决策变量x_ij:若第i个人选择第j项任务则x_ij=1;否则，x_ij=0；

模型为：

model:
sets:
task/1..4/:t;
man/1..4/:m;
link(man,task):c,x;
endsets
data:
c=2 15 13 4
10 4 14 15
9 14 16 13
7 8 11 9;
enddata
[obj]min=@sum(link:c*x);
@for(task(j):@sum(man(i):x(i,j))=1);
@for(man(i):@sum(task(j):x(i,j))=1);
@for(link:@bin(x));
end
求得：最优指派为：m1--t4,m2--t2,m3--t1,m4--t3

最优值为:28。

 

六、装配线平衡模型（0-1规划问题）

11 件任务（A—K）分配到 4 个工作站（1—4），任务的优先次序如下图，每件任务所花费的时间如下表。目标是为每个工作站分配加工任务，尽可能使每个工作站执行相同的任务量，其最终装配线周期最短。

任务 A  B  C  D  E  F  G  H  I  J  K

时间 45 11 9 50 15 12 12 12 12  8  9 

T(i):为完成第i项任务需要的时间。

 

SETS:

TASK/ A B C D E F G H I J K/: T; !任务集合，有一个完成时间属性 T;

PRED( TASK, TASK)/ A,B B,C C,F C,G F,J G,J

J,K D,E E,H E,I H,J I,J /; !任务之间的优先关系集合（A 必须完成才能开始 B，等等）;

STATION/1..4/; ! 工作站集合;

TXS( TASK, STATION): X;! X 是派生集合 TXS 的一个属性。如果 X（I，K）＝1，则表

示第 I 个任务指派给第 K 个工作站完成;

ENDSETS

DATA:

T = 45 11 9 50 15 12 12 12 12 8 9; !任务 A B C D E F G H I J K 的完成时间;

ENDDATA

@FOR( TASK( I): @SUM( STATION( K): X( I, K)) = 1); !每一个作业必须指派到一个工

作站;

@FOR( PRED( I, J): @SUM( STATION( K):  X(I, K))-@SUM( STATION( K):  X(J, 

K) )>=0) ; !对于每一个存在优先关系的作业对(I,J)来说，I先J后安排;

@FOR( STATION( K):@SUM( TXS( I, K): T( I) * X( I, K)) <= CYCTIME); !对于每一个

工作站来说，其花费时间必须不大于装配线周期;

MIN = CYCTIME; !目标函数是最小化转配线周期;

@FOR( TXS: @BIN( X)); !指定 X(I,J) 为 0/1 变量;

END

 

解得最短周期为50.

分配情况为：A-1，B-3，C-4，D-2，E-3，F-4，G-4，H-3，I-3，J-4，K-4.

 

七、选址问题

某海岛上有12个主要的居民点，每个居民点的位置（用平面坐标x,y表示，距离单位：km）和居住的人数（r）如下表所示。现在准备在海岛上建一个服务中心为居民提供各种服务，那么服务中心应该建在何处？
x 0 8.20 0.50 5.70 0.77 2.87 4.43 2.58 0.72 9.76 3.19 5.55

y 0 0.50 4.90 5.00 6.49 8.76 3.26 9.32 9.96 3.16 7.20 7.88

r 600 1000 800 1400 1200 700 600 800 1000 1200 1000 1100

设建在（a，b）处最合理。建立模型：

MODEL:
SETS:
VAR/1..12/:X,Y,R;
ENDSETS
DATA:
X=0 8.20 0.50 5.70 0.77 2.87 4.43 2.58 0.72 9.76 3.19 5.55;
Y=0 0.50 4.90 5.00 6.49 8.76 3.26 9.32 9.96 3.16 7.20 7.88;
R=600 1000 800 1400 1200 700 600 800 1000 1200 1000 1100;
ENDDATA
MIN=@SUM(VAR:@SQRT((X-A)^2+(Y-B)^2)*R);
END

求得:(a,b)=(3.601028,6.514223),最小值为：44236.04。

 

八、婚配问题：

10对年龄相当的青年，任意一对男女青年配对的概率pij见下表。试给出一个配对方案，使总的配对概率最大。

    w1        w2        w3        w4        w5        w6        w7        w8        w9        w10
m1 0.5828    0.2091    0.4154    0.2140    0.6833    0.4514    0.6085    0.0841    0.1210  0.2319
m2 0.4235    0.3798    0.3050    0.6435    0.2126    0.0439    0.0158    0.4544    0.4508  0.2393
m3 0.5155    0.7833    0.8744    0.3200    0.8392    0.0272    0.0164    0.4418    0.7159  0.0498
m4 0.3340    0.6808    0.0150    0.9601    0.6288    0.3127    0.1901    0.3533    0.8928  0.0784
m5 0.4329    0.4611    0.7680    0.7266    0.1338    0.0129    0.5869    0.1536    0.2731  0.6408
m6 0.2259    0.5678    0.9708    0.4120    0.2071    0.3840    0.0576    0.6756    0.2548  0.1909
m7 0.5798    0.7942    0.9901    0.7446    0.6072    0.6831    0.3676    0.6992    0.8656  0.8439
m8 0.7604    0.0592    0.7889    0.2679    0.6299    0.0928    0.6315    0.7275    0.2324  0.1739
m9 0.5298    0.6029    0.4387    0.4399    0.3705    0.0353    0.7176    0.4784    0.8049  0.1708

m10 0.6405    0.0503    0.4983    0.9334    0.5751    0.6124    0.6927    0.5548    0.9084  0.9943

取xx_ij为0-1型决策变量。

模型为：

model:
sets:
man/m1..m10/;
woman/w1..w10/;
link(man,woman):p,x;
endsets
data:
p=0.5828    0.2091    0.4154    0.2140    0.6833    0.4514    0.6085    0.0841    0.1210    0.2319
  0.4235    0.3798    0.3050    0.6435    0.2126    0.0439    0.0158    0.4544    0.4508    0.2393
  0.5155    0.7833    0.8744    0.3200    0.8392    0.0272    0.0164    0.4418    0.7159    0.0498
  0.3340    0.6808    0.0150    0.9601    0.6288    0.3127    0.1901    0.3533    0.8928    0.0784
  0.4329    0.4611    0.7680    0.7266    0.1338    0.0129    0.5869    0.1536    0.2731    0.6408
  0.2259    0.5678    0.9708    0.4120    0.2071    0.3840    0.0576    0.6756    0.2548    0.1909
  0.5798    0.7942    0.9901    0.7446    0.6072    0.6831    0.3676    0.6992    0.8656    0.8439
  0.7604    0.0592    0.7889    0.2679    0.6299    0.0928    0.6315    0.7275    0.2324    0.1739
  0.5298    0.6029    0.4387    0.4399    0.3705    0.0353    0.7176    0.4784    0.8049    0.1708
  0.6405    0.0503    0.4983    0.9334    0.5751    0.6124    0.6927    0.5548    0.9084   0.9943;
enddata
max=@prod(man(i):@sum(woman(j):p(i,j)*x(i,j)));
@for(woman(j):@sum(link(i,j):x(i,j))=1);
@for(man(i):@sum(link(i,j):x(i,j))=1);
@for(link:@bin(x));
end

求解结果：m1-w5,m2-w8,m3-w2,m4-w4,m5-w7,m6-w3,m7-w6,m8-w1,m9-w9,m10-w10.

最大值为0.055.

 

九、护士值班安排问题

    某医院，从周一到周日都要有人值班，每天至少需要的护士如表。要求每个护士每周连续上班5天，试问该医院至少需要多少护士？并给出上班安排计划。

周   1    2    3    4    5     6     7

人   20   16   13   16   19    14    12

 

取决策变量star（i）:周i开始值班的人数；

目标函数:min  sum[star(i)]（i=1,2,3，...,7)

约束条件：连续工作五天，周j值班的人数>=required（j）（j=1,2,3,...,7)

 

model:

sets:

days/mon..sun/: required,start;

endsets

data:

!每天所需的最少职员数;

required = 20 16 13 16 19 14 12;

enddata

!最小化每周所需职员数;

min=@sum(days: start);

@for(days(J):@sum(days(I) | I #le# 5:start(@wrap(J+2+I,7))) >= required(J));

end

解得： 总共需要22人，周一8人开始值班，周二2人，周三0人，周四6人，周五3人，周六3人，周日0人。

 

十、填数问题

     分别把1,2,…,16填到图示的16个圈内，使得每个三角形上的所有圈内的数的和为81（共4个三角形）。
决策变量:e_ij=1,第i个圈填数a_j;e_ij=0,第i个圈不填数a_j。a_j=j,j=1,2,3,...,16。

模型：

model:

sets:

number/1..16/:a;

link(number,number):e;

tri1(number)/1 2 3 4 5 6 7 8 9/;

tri2(number)/1 2 3 4 16 15 12 11 10/;

tri3(number)/4 5 6 7 14 13 12 15 16/;

tri4(number)/7 8 9 1 10 11 12 13 14/;

endsets

data:

a=1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16;

enddata

[obj]max=@sum(link(i,j):e(i,j)*a(j));

@for(number(i):@sum(link(i,j):e(i,j))=1);

@for(number(j):@sum(link(i,j):e(i,j))=1);

@for(link(i,j):@bin(e(i,j)));

@sum(number(j):@sum(tri1(i):e(i,j)*a(j)))=81;

@sum(number(j):@sum(tri2(i):e(i,j)*a(j)))=81;

@sum(number(j):@sum(tri3(i):e(i,j)*a(j)))=81;

@sum(number(j):@sum(tri4(i):e(i,j)*a(j)))=81;

@sum(link(i,j):e(i,j)*a(j))=136;

end

红色的那句程序可以去掉，也可以为：min=@sum(link(i,j):e(i,j)*a(j))，但求的结果不同，结果都符合要求。

编号1~16的圆圈的填数结果至少有3种：

（1）12 11 1 10 7 8 14 13 5 9 4 16 2 6 15 3

（2）14 3 5 15 8 7 13 4 12 6 11 10 9 2 16 1

（3）14 11 4 15 9 8 13 2 5 16 3 10 12 6 1 7

为了求得更多的解，可以约束编号1~16的圆圈的填数结果不为以上3种结果。

model:

sets:

number/1..16/:a;

link(number,number):e;

tri1(number)/1 2 3 4 5 6 7 8 9/;

tri2(number)/1 2 3 4 16 15 12 11 10/;

tri3(number)/4 5 6 7 14 13 12 15 16/;

tri4(number)/7 8 9 1 10 11 12 13 14/;

yueshu1:c1;

yueshu2:c2;

yueshu3:c3;

endsets

data:

a=1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16;

c1=12 11 1 10 7 8 14 13 5 9 4 16 2 6 15 3;

c2=14 3 5 15 8 7 13 4 12 6 11 10 9 2 16 1;

c3=14 11 4 15 9 8 13 2 5 16 3 10 12 6 1 7;

enddata

[obj]min=@sum(number(i):@sum(number(j):e(i,j)*a(j)));

@for(number(i):@sum(number(j):e(i,j))=1);

@for(number(j):@sum(link(i,j):e(i,j))=1);

@for(link(i,j):@bin(e(i,j)));

@sum(number(j):@sum(tri1(i):e(i,j)*a(j)))=81;

@sum(number(j):@sum(tri2(i):e(i,j)*a(j)))=81;

@sum(number(j):@sum(tri3(i):e(i,j)*a(j)))=81;

@sum(number(j):@sum(tri4(i):e(i,j)*a(j)))=81;

@sum(link(i,j):e(i,j)*a(j))=136;

@sum(yueshu1(j):@sum(link(i,j):e(i,j)))<16;

@sum(yueshu2(j):@sum(link(i,j):e(i,j)))<16;

@sum(yueshu3(j):@sum(link(i,j):e(i,j)))<16;

end

 

解得：（4）12 15 11 10 3 8 16 1 5 4 7 14 9 13 2 6

可以继续下去：······

 

十一、数独求解

在琳琅在线数独随机截取一至尊级题目如下：

用Lingo求解这道题。

Lingo代码：

model:
sets:
row/1..9/;
col/1..9/;
num/1..9/;
var(row,col,num):x;
endsets
@for(col(j):@for(num(k):@sum(var(i,j,k):x)=1));!约束每一列填的数字不重复;
@for(row(i):@for(num(k):@sum(var(i,j,k):x)=1));!约束每一行填的数字不重复;
@for(row(i):@for(col(j):@sum(var(i,j,k):x)=1));!约束每一个格子必须填入数字;
!约束9个块里的数字不重复;
@for(num(k):@sum(row(i)|i#ge#1#and#i#le#3:
@sum(col(j)|j#ge#1#and#j#le#3:x(i,j,k)))=1;

@sum(row(i)|i#ge#1#and#i#le#3:
@sum(col(j)|j#ge#4#and#j#le#6:x(i,j,k)))=1;

@sum(row(i)|i#ge#1#and#i#le#3:
@sum(col(j)|j#ge#7#and#j#le#9:x(i,j,k)))=1;

@sum(row(i)|i#ge#4#and#i#le#6:
@sum(col(j)|j#ge#1#and#j#le#3:x(i,j,k)))=1;

@sum(row(i)|i#ge#4#and#i#le#6:
@sum(col(j)|j#ge#4#and#j#le#6:x(i,j,k)))=1;

@sum(row(i)|i#ge#4#and#i#le#6:
@sum(col(j)|j#ge#7#and#j#le#9:x(i,j,k)))=1;

@sum(row(i)|i#ge#7#and#i#le#9:
@sum(col(j)|j#ge#1#and#j#le#3:x(i,j,k)))=1;

@sum(row(i)|i#ge#7#and#i#le#9:
@sum(col(j)|j#ge#4#and#j#le#6:x(i,j,k)))=1;

@sum(row(i)|i#ge#7#and#i#le#9:
@sum(col(j)|j#ge#7#and#j#le#9:x(i,j,k)))=1);
!已知的条件;
x(1,6,7)=1;x(1,7,4)=1;
x(2,1,3)=1;x(2,4,8)=1;x(2,9,1)=1;
x(3,5,6)=1;x(3,8,9)=1;
x(4,5,4)=1;x(4,8,7)=1;
x(5,3,5)=1;x(5,4,1)=1;x(5,9,2)=1;
x(6,6,9)=1;x(6,7,6)=1;
x(7,1,2)=1;x(7,3,3)=1;
x(8,2,6)=1;x(8,3,8)=1;x(8,4,5)=1;
x(9,1,5)=1;x(9,2,1)=1;x(9,9,3)=1;
@for(var:@bin(x));!0,1约束;
end

结果：