博客中所有内容均来源于自己学习过程中积累的经验以及对yalmip官方文档的翻译:https://yalmip.github.io/tutorials/
之前的博客简单介绍了约束条件的定义方法,接下来将对其进行详细介绍。
首先简单复习一下:
1.定义约束条件可以使用矩阵拼接,或者’+’号,两种方式是等价的;
2.Yalmip中的约束可以写成连续的不等式形式;
3.Yalmip中不支持严格的不等号(<或>),一旦出现,就会报错,然后弹出“猫猫图”。
下面将继续讲解约束条件定义的方法与技巧。
1.约束的普通形式与矩阵形式
例1:请使用Yalmip工具箱写出下列优化问题的约束条件
采用常规的思维,可以分别写出每个约束条件,然后采用拼接的方式将所有约束条件组合起来,matlab代码如下:
x = sdpvar(3,1);C = [sum(x) <= 12 ,...5*x(1) + x(3) <= 9 ,...3*x(1) + 2*x(2) >= 6 ,...x(2) + 2*x(3) >= 1 ,...x >= 0];
但是在求解优化问题的过程中,很多时候需要我们写出优化问题的对偶形式或KKT条件,这时候如果可以写成紧凑的矩阵形式,转换将会更加容易。例1的约束条件写成紧凑的矩阵形式为:
将约束条件写作矩阵形式时,代码如下:
x = sdpvar(3,1);
A = [-1 0 0;0 -1 0;0 0 -1;1 1 1;5 0 1;-3 -2 0;0 -1 -2];
b = [0;0;0;12;9;-6;-1];
C = [A*x <= b];
两种写法是完全等价的,下面的代码分别采用两种不同的约束定义方式对优化问题进行求解,可以得到相同的优化结果:
clc
clearx = sdpvar(3,1);
A = [-1 0 0;0 -1 0;0 0 -1;1 1 1;5 0 1;-3 -2 0;0 -1 -2];
b = [0;0;0;12;9;-6;-1];
C1 = [A*x <= b];
C2 = [sum(x) <= 12 ,...5*x(1) + x(3) <= 9 ,...3*x(1) + 2*x(2) >= 6 ,...x(2) + 2*x(3) >= 1 ,...x >= 0];
objective = [1 2 3]*x;
optimize(C1 , objective);
objective1 = value(objective)
x1 = value(x)
optimize(C2 , objective);
objective2 = value(objective)
x2 = value(x)
运行结果:
objective1 =3.3333x1 =1.33331.00000objective2 =3.3333x2 =1.33331.00000
在实际编程过程中,两种形式各有优劣,但完全等价,可以根据自己的需要采用一般形式或矩阵形式定义约束条件,只要自己感觉更方便即可。
2.小技巧
2.1查看约束变量的参数
在定义约束条件之后,在工作区双击对应的约束变量,或者在命令行输入约束变量再回车,可以查看约束的列表,并检查约束是否已有效定义了,例如,例1中用两种不同方式定义的约束条件列表分别如下:
在该列表中,ID一列表示约束的编号,每个约束都是用’,’或者’+’隔开的,因此约束变量C1中仅有1个约束,约束变量C2中共有5条约束。Constraint一列表示约束的类型和数目,Element-wise inequality表示该约束条件为不等式约束,如果是等式约束,则显示‘Equality constraint’,其他一些常见的约束类型如表1所示。约束类型后面显示的是约束的数目,例如约束变量C1共包含七条约束,因此显示为7×1的约束,而约束变量C2中的第五条约束x>=0,实际上包含x1>=0,x2>=0,x3>=0共三条约束,因此显示为3×1的约束。
表1 常见的约束类型
约束类型 | 含义 |
Matrix inequality | 矩阵不等式约束 |
Second order cone constraint | 二阶锥约束 |
Rotated Lorentz constraint | 旋转锥约束 |
Binary constraint | 二进制约束 |
Eigenvalue constraint | 特征值约束 |
KYP constraint | Kalman Yakubovich Popov 约束 |
Sum-of-square constraint | 平方和约束 |
Logic constraint | 逻辑约束 |
Semi-continuous variable | 半连续的变量约束 |
Semi-integer variable | 半整数的变量约束 |
Vectorized second-order cone constraints | 矢量化二阶锥约束 |
2.2尽量避免使用循环语句
在Yalmip工具箱中,大部分matlab内置函数都可以用于sdpvar变量。之前我写过几篇博客说明如何在matlab代码中尽量避免循环的使用,以提升代码运行效率,因此也可将这种思想用于约束条件的定义过程中。下面以Yalmip工具箱官方文档中给出的数独求解的算例(https://yalmip.github.io/example/sudoku/)进行说明。
例2:使用Yalmip工具箱求解图1所示的数独问题
图1 一个数独问题
为了求解这个问题,可以定义一个9×9×9的三维0-1变量A(i,j,k),当A(i,j,k)取值为1时,表示该数独的第i行第k列的取值为k,为了求解上述数独问题,需要写成初始的数独矩阵,并定义决策变量A:
A0 = [
0 4 7 0 5 0 0 0 8;
6 0 5 0 3 0 2 0 1;
0 0 0 7 0 6 0 3 0;
0 0 6 0 7 0 0 2 4;
9 0 0 8 0 4 0 0 6;
4 5 0 0 1 0 9 0 0;
0 1 0 5 0 2 0 0 0;
2 0 8 0 4 0 5 0 3;
5 0 0 0 9 0 7 1 0;
];A = binvar(9,9,9,'full');
首先,我们需要把变量A中的元素取值和初始矩阵A0对应起来,代码如下:
for i = 1:9for j = 1:9if S(i,j)F = [F, A(i,j,A0(i,j)) == 1];endend
end
这也是常规使用循环语句的思维,但是,如果使用find和sub2ind命令,可以避免循环语句的使用。find函数比较常用,也比较简单,这里不再介绍,重点介绍一下sub2ind。该函数用于将下标向量组转为线性索引,具体的含义我们可以通过一个例子来理解。假设在matlab中有一个矩阵3×3的矩阵S,我们如果想取出矩阵S第3行第2列的元素,可以采用下标的形式写成S(3,2),如果采用线性索引的话,也可以写成A(6)。下标和线性索引的关系如图2所示:
图2 matlab二维数组中下标和线性索引的关系
sub2ind就可以把一组下标转为一组线性索引,以避免循环的使用,标准语法为:
ind = sub2ind(sz,row,col)
例如,我们要取出矩阵S的元素S(1,2),S(1,3),S(2,2)和S(3,2),也就是要求出线性索引向量[4 5 6 7],采用sub2ind函数即可解决:
row = [1 2 3 1];
col = [2 2 2 3];
sz = [3 3];
ind = sub2ind(sz,row,col)
输出结果为:
ind =4 5 6 7
因此,使用find与sub2ind,可以将上述代码改写成:
[i,j,k] = find(A0);
F = [F, A(sub2ind(size(A),i,j,k)) == 1];
这样做就可以避免使用双层循环。另外,数独游戏的规则是,每一行、每一列以及每一个3×3的子块的数字不能重复,其中每行每列数字不能重复的约束比较简单,代码如下:
F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];
要保证每个3×3的子块数字不能重复,可以采用循环的方式,代码如下:
for m = 1:3for n = 1:3for k = 1:9s = sum(sum(A((m-1)*3+(1:3),(n-1)*3+(1:3),k))); F = [F, s == 1];endend
end
为了尽量避免使用循环语句,可以将代码改写如下:
for k = 1:9F = [F, sum(sum(A(i1:i1+2 , j1:j1+2 , :))) == 1];
end
这样写可以只用一层循环,具体原理比较简单,这里不再过多讲解。另外,对于数独问题,我们只需要求出可行解,并没有优化目标,因此使用optimize函数时可以只输入约束条件。
sol = optimize(F);
求解成功后,可以使用一个矩阵Z来存储最终的数独求解结果:
Z = 0;
for i = 1:9Z = Z + i*value(A(:,:,i));
end
Z
为了避免循环语句,也可以写成:
Z = value(reshape(A,9,81)*kron((1:9)',eye(9)))
其中,reshape函数用于改变矩阵的形状,并返回一个新的矩阵。在这条语句中,reshape(A,9,81)表示将矩阵A重塑为一个9行81列的矩阵。kron函数用于计算两个矩阵的张量积,在这条语句中,kron((1:9)',eye(9))可以得到一个81行9列的矩阵,其中对应位置的元素由(1:9)'和单位矩阵的元素相乘得到。将矩阵reshape(A,9,81)与矩阵kron((1:9)',eye(9)),得到的就是9×9的矩阵,也就是数独问题的解。
综上所示,使用常规思维+更多的循环语句求解数独问题的完整代码如下:
方法1:放肆使用循环语句
clc
cleartic
A0 = [
0 4 7 0 5 0 0 0 8;
6 0 5 0 3 0 2 0 1;
0 0 0 7 0 6 0 3 0;
0 0 6 0 7 0 0 2 4;
9 0 0 8 0 4 0 0 6;
4 5 0 0 1 0 9 0 0;
0 1 0 5 0 2 0 0 0;
2 0 8 0 4 0 5 0 3;
5 0 0 0 9 0 7 1 0;
];A = binvar(9,9,9,'full');
F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];for i = 1:9for j = 1:9if A0(i,j)F = [F, A(i,j,A0(i,j)) == 1];endend
endfor m = 1:3for n = 1:3for k = 1:9s = sum(sum(A((m-1)*3+(1:3),(n-1)*3+(1:3),k))); F = [F, s == 1];endend
enddiagnostics = optimize(F);Z = 0;
for i = 1:9Z = Z + i*value(A(:,:,i));
end
Z
toc
方法2:尽量减少循环语句
clc
cleartic
A0 = [
5 3 0 0 7 0 0 0 0;
6 0 0 1 9 5 0 0 0;
0 9 8 0 0 0 0 6 0;
8 0 0 0 6 0 0 0 3;
4 0 0 8 0 3 0 0 1;
7 0 0 0 2 0 0 0 6;
0 6 0 0 0 0 2 8 0;
0 0 0 4 1 9 0 0 5;
0 0 0 0 8 0 0 7 9;
];A = binvar(9,9,9,'full');F = [sum(A,1) == 1, sum(A,2) == 1, sum(A,3) == 1];i1 = [1,1,1,4,4,4,7,7,7];
j1 = [1,4,7,1,4,7,1,4,7];
for k = 1:9F = [F, sum(sum(A(i1:i1+2 , j1:j1+2 , :))) == 1];
end[i,j,k] = find(A0);
F = [F, A(sub2ind(size(A),i,j,k)) == 1];diagnostics = optimize(F);Z = value(reshape(A,9,81)*kron((1:9)',eye(9)))
toc
两种方式的运行结果完全一致:
Z =3 4 7 2 5 1 6 9 86 8 5 4 3 9 2 7 11 2 9 7 8 6 4 3 58 3 6 9 7 5 1 2 49 7 1 8 2 4 3 5 64 5 2 6 1 3 9 8 77 1 3 5 6 2 8 4 92 9 8 1 4 7 5 6 35 6 4 3 9 8 7 1 2
但在运行时间上,两种方法存在差异,某次测试中,方法1的运行时间为0.638127秒,方法2的运行时间为0.473180秒,避免循环的使用提高了25.85%的运行效率。
减少循环语句的使用,在小规模问题中效率提升可能不是很明显,但如果优化问题比较复杂,约束也很复杂,如果可以把尽量减少循环语句的思想应用到实际编程中,效率提升可能是非常大的。虽然思考如何减少循环语句可能要废很多脑细胞,但我觉得应该会挺值得。(我自己经历过一个优化问题求解需要8个小时左右,每次等待结果、调试代码都非常痛苦,将近半个月都没有调好,后来想通了,花了两三天的时间优化代码,把运行时间缩短到了1个半小时左右,很快就调试成功了)
3.测试题
3.1测试1
使用矩阵形式的约束条件求解下列优化问题:
max z=3x1+x2
s.t. x1-x2≥-2
x1-2x2≤3
3x1+2x2≤14
x1≥0, x2≥0
3.2测试2
使用Yalmip工具箱求下列数独问题的解:
3.3测试3
改写下面的代码以避免使用循环语句,并对比两种方式的代码运行时间。
clc
clear
yalmip('clear')tic
n = 50;
a = rand(n,n,n);
x = sdpvar(n,n,n);
z = 0;
for k = 1:nfor kk = 1:nfor kkk = 1:nz = z + a(k,kk,kkk) * x(k,kk,kkk);endend
end
toc
3.4测试题参考答案
第三章测试题的参考答案可以从下面的链接中获取:
https://download.csdn.net/download/weixin_44209907/88062937