赞
踩
目录
1.2方阵默认为对称阵。若需要定义一个完全参数化(即不一定对称)的方阵,还需要输入第三个参数:
1.3如果在上述代码中省略句尾的分号,或直接在命令行输入“P”后回车,可以查看此矩阵的性质。结果如下所示:
yalmip是由Lofberg开发的一种免费的优化求解工具,其最大特色在于集成许多外部的最优化求解器,形成一种统一的建模求解语言,提供了Matlab的调用API,减少学习者学习成本。
首先,yalmip是一个matlab的工具包,通过matlab实现各种操作和调用。
其次,它是一个建模工具,甚至可以称为一种“语言”,通过这种“语言”来描述模型,然后再调用其他求解器(如gurobi、cplex等)来求解模型。相当于一个将“yalmip语言”转换成其他求解器“语言”的语言转换器。
不同的求解器有不同的专用语言,学习多个语言即冗余又浪费精力,所以,yalmip的珍贵之处就体现出来了。
更为可贵的是,yalmip真正实现了建模和算法二者的分离,它提供了一种统一的、简单的建模语言,针对所有的规划问题,都可以用这种统一的方式建模;至于用哪种求解算法,你只需要通过一次简单的参数配置指定就可以了,甚至不用你指定,yalmip会自动为你选择最适合的算法。
有了yalmip,你不再需要针对每一种工具包去学习特定的建模语言(比如用cplex要专门学习cplex的建模语言,用lingo要专门学习lingo的建模语言,还有GLPK、lpsolve、Matlab自带的求解器等等,如果每一种求解器都要学习新的建模语言的话,这个工作量是可想而知的)。相反,如果你选择使用yalmip,那么你只需要学习yalmip一种建模语法,因为yalmip真正实现了建模和算法的分离,所有的问题都可以用统一的方法建模,如果需要使用不同的求解器,只需要一句简单的配置即可。因此,yalmip不仅仅是一个线性规划求解器,更强大的地方在于,它提供了一个统一的建模平台,支持现有的几乎所有的求解算法。有了yalmip,一切都变得简单起来。
以上摘自博客《yalmip + lpsolve + matlab 求解混合整数线性规划问题(MIP/MILP)》
这里以MATLAB的安装方式为例
官网下载https://yalmip.github.io/ download 即可 ,可得到YALMIP-master .zip
解压至matlab/toolbox
matlab中设置路径,注意:要将压缩包内的子文件夹都加入路径,选择“添加并包含子文件夹”来添加路径;
检查是否配置成功:matlab中调用yalmitest命令,查看所有支持的求解器已经他们的安装状态。
最后键入which sdpvar
命令,显示sdpvar
路径则安装成功
以下参考博文 https://www.jianshu.com/p/e1c45b3d8d8a
yalmip一共有三种方式创建决策变量,分别为:
- sdpvar-创建实数型决策变量
- intvar-创建整数型决策变量
- binvar-创建0/1型决策变量
不过值得注意的是,在创建n*n的决策变量时,yalmip默认是对称方阵。
所以要创建非对称方针时,需要这样写:xxxvar(n,n,'full')
比起matlab自带的各种优化函数所要写明的约束条件,yalmip的约束条件写起来是非常舒适直观的。符合人类直觉
比如要写入0<=x1+x2+x3<=1。
那么可以这样写:
- % 创建决策变量
- x = sdpvar(1,3);
- % 添加约束条件
- C = [0<=x(1)+x(2)+x(3)<=1];
关于参数设置,我们大多数是用来设置求解器solver的,当然还有其它的选项,可以通过doc sdpsettings查看。
最后就是求解问题了。
首先要明确求解目标z,yalmip默认是求解最小值问题,所以遇到求解最大值的问题,只需要在原问题的基础上添加一个负号即可。
求解调用格式:
- 求解调用格式:
-
- optimize(target,constraints,opstions)
- clear;clc;close all
- x = sdpvar(1,3);
- z = 2*x(1) + 3*x(2) + x(3);
- c = [x(1) + 4*x(2) + 2*x(3) >= 8
- 3*x(1) + 2*x(2) >= 6
- x(1), x(2), x(3) > 0];
-
- result = optimize(c,z);
-
- if result.problem == 0 % 求解成功
- xresult = value(x)
- zresult = value(z)
- else
- disp('求解出错')
- end
运行结果如下
- CPXPARAM_MIP_Display 1
- Tried aggregator 1 time.
- LP Presolve eliminated 3 rows and 0 columns.
- Reduced LP has 2 rows, 3 columns, and 5 nonzeros.
- Presolve time = 0.02 sec. (0.00 ticks)
-
- Iteration log . . .
- Iteration: 1 Dual objective = 4.000000
-
- xresult =
-
- 2 0 3
-
-
- zresult =
-
- 7
当然该问题也可以直接用 matlab自带的linprog进行求解,代码如下
- c = [2 3 1]
- a = [1 4 2;
- 3 2 0];
- b = [8;6];
- [x,z] = linprog(c,-a,-b,[],[],zeros(3,1))
结果如下,与上面的一致。
- Optimal solution found.
-
-
- x =
-
- 2.0000
- 0
- 3.0000
-
-
- z =
-
- 7.0000
- % 清除工作区
- clear;clc;close all;
- % 创建决策变量
- x = sdpvar(1,2);
- % 添加约束条件
- C = [
- x(1) + x(2) >= 2
- x(2)-x(1) <=1
- x(1)<=1
- ];
- % 配置
- ops = sdpsettings('verbose',0,'solver','lpsolve');
- % 目标函数
- z = -(x(1)+2*x(2))/(2*x(1)+x(2)); % 注意这是求解最大值
- % 求解
- reuslt = optimize(C,z);
- if reuslt.problem == 0 % problem =0 代表求解成功
- value(x)
- -value(z) % 反转
- else
- disp('求解出错');
- end
求解结果:
- %求解如下非线性规划问题
-
-
- % 清除工作区
- clear;clc;close all;
- % 创建决策变量
- x = sdpvar(1,2);
- % 添加约束条件
- C = [
- x(1) + x(2) >= 2
- x(2)-x(1) <=1
- x(1)<=1
- ];
- % 配置
- ops = sdpsettings('verbose',0,'solver','cplex');
- % 目标函数
- z = -(x(1)+2*x(2))/(2*x(1)+x(2)); % 注意这是求解最大值
- % 求解
- reuslt = optimize(C,z);
- if reuslt.problem == 0 % problem =0 代表求解成功
- value(x)
- -value(z) % 反转
- else
- disp('求解出错');
- end
求解结果:
- clear;clc;close all;
- %定义变量
- x=sdpvar(2,1);
-
- %约束条件
- constraint=[];
- constraint=[constraint,x(1)+x(2)>350];
- constraint=[constraint,x(1)>100];
- constraint=[constraint,2*x(1)+x(2)<600];
- constraint=[constraint,x(2)>0];
- %求解
- %ops = sdpsettings('solver','cplex','verbose',1);
- % ops.cplex.display='on';
- % ops.cplex.timelimit=600;
- % ops.cplex.mip.tolerances.mipgap=0.001;
-
-
-
- % 诊断求解可行性
- % disp('开始求解')
- % diagnostics=optimize(constraint,obj,ops);
- % if diagnostics.problem==0
- % disp('Solver thinks it is feasible')
- % elseif diagnostics.problem == 1
- % disp('Solver thinks it is infeasible')
- % pause();
- % else
- % disp('Timeout, Display the current optimal solution')
- % end
-
-
- % 配置
- ops = sdpsettings('solver','cplex'); %配置求解方法为调用 CPLEX
-
-
- %目标函数
- obj=2*x(1)+3*x(2);
-
- % 求解
- reuslt = optimize(constraint,obj,ops); %Yalmip优化求解的命令
- if reuslt.problem == 0 % problem =0 代表求解成功
- value(x)
- value(obj)
- else
- disp('求解出错');
- end
求解结果:
在YALMIP中最重要的命令是sdpvar,用于定义决策变量。
P = sdpvar(n,m);
P = sdpvar(3,3,'full');
Linear matrix variable 3x3 (full, real, 9 variables)
第三个参数也可以用于获得一系列预定义的变量,如Toeplitz、Hankel、对角、对称和斜对称矩阵等(更多细节请参考sdpvar)。当然,也可以用常规方法来定义向量,如下所示:
- x = sdpvar(n,1); % 向量
- D = diag(x) ; % 对角矩阵
- H = hankel(x); % Hankel矩阵
- T = toeplitz(x); % Toeplitz矩阵
标量可以通过如下三种方法来定义:
- x = sdpvar(1,1); y = sdpvar(1,1);
- x = sdpvar(1); y = sdpvar(1);
- sdpvar x y
sdpvar对象在MATLAB中的使用方法和任何其他变量一样,大多数函数都是重载的。因此,以下命令是有效的:
- P = sdpvar(3,3) + diag(sdpvar(3,1));
- X = [P P;P eye(length(P))] + 2*trace(P);
- Y = X + sum(sum(P*rand(length(P)))) + P(end,end)+hankel(X(:,1));
-
- % eye()函数用于返回单位矩阵,如eye(n)将返回n*n单位矩阵
- % trace()函数用于求矩阵的迹
- % end函数返回下标的最大值
在码代码时,如果对YALMIP不太熟悉,建议养成查看表达式和变量的习惯,看看它们的属性是否与期望相符。举个栗子,若定义正确,上文代码所得的X应该是一个6x6实对称矩阵。根据下文对X的查看结果可得,它确实是一个由9个变量构成的对称6x6矩阵。
- >> X
- Linear matrix variable 6x6 (symmetric, real, 9 variables)
在某些情况下,可以使用多维变量来简化编码。YALMIP支持两种不同的构造:元胞数组和多维sdpvar对象。元胞数组格式只是以下代码的抽象化:
- for i = 1:5
- X{i} = sdpvar(2,3);
- end
通过在sdpvar中设置向量维数,上文的元胞数组还可以用下述方法设置:
X = sdpvar([2 2 2 2 2],[3 3 3 3 3]);
此元胞数组可以在MATLAB中正常使用。但这种定义方法的缺陷在于变量X不能直接用作标准的sdpvar对象(在MATLAB中,plus等运算符在单元格上不可重载)。因此,可以使用完全通用的多维sdpvar创建一个等价对象:
X = sdpvar(2,3,5);
通过这种方法,我们可以用标准的MATLAB代码直接对这个对象进行操作:
- Y = sum(X,3)
- X(:,:,2)
根据标准YALMIP语法,如果前两个维度相同,那么前两个面就是对称的。要创建完全参数化的高维变量,需要给定一个“full”标志:
X = sdpvar(2,2,2,2,'full');
有关多维变量的示例请查看Sudoku example。
我们可以通过“创建+连接”的方式定义约束集合。约束的含义取决于上下文:如果左右两侧都是Hermitian矩阵,则约束将按照正定性进行解释,否则将逐元素进行解释。因此,定义一个对称矩阵和正定性约束的方法如下:
- n = 3;
- P = sdpvar(n,n);
- C = [P>=0];
定义具有正元素的对称矩阵示例如下:
- P = sdpvar(n,n);
- C = [P(:)>=0];
注意,这样就两次定义了非对角约束。一个好的半定规划(Semi-definite Programming,SDP)求解器可能会在预处理过程中检测到这一点并简化模型,不过我们也可以使用常规方法来手动定义特殊元素,如下所示:
C = [triu(P)>=0]; % triu()函数用于抽取上三角矩阵
或者:
C = [P(find(triu(ones(n))))>=0];
根据上述规则,可以使用>=运算符定义具有正元素的非平方矩阵(即非对称矩阵),如下所示:
- P = sdpvar(n,2*n);
- C = [P>=0];
同样的方法也可以用于定义具有正元素的完全参数化的矩阵:
- P = sdpvar(n,n,'full');
- C = [P>=0];
一系列约束可以通过“添加”或“连接”来定义,如下所示:
- P = sdpvar(n,n);
- C = [P>=0] + [P(1,1)>=2];
- C = [P>=0, P(1,1)>=2];
定义约束时所涉及的表达式中可以包括任意sdpvar对象,等式约束(==)和不等式约束(<=)也都可以使用,如下所示:
C = [P>=0, P(1,1)<=2, sum(sum(P))==10
在工作区双击约束,或在命令行输入约束名后回车,可以查看约束列表,并检查约束是否已有效定义了。正如前文关于查看表达式和变量的建议一样,推荐大家及时查看约束,以保证它们的属性与期望相符。查看上述约束所得结果示例如下:
- ++++++++++++++++++++++++++++++++++++++
- | ID| Constraint|
- ++++++++++++++++++++++++++++++++++++++
- | #1| Matrix inequality 3x3|
- | #2| Element-wise inequality 1x1|
- | #3| Equality constraint 1x1|
- ++++++++++++++++++++++++++++++++++++++
以下代码可以实现多个约束的同时定义:
F = [0 <= P(1,1) <= 2];
严格不等式是不能使用的(将收到YALMIP发出的警告),因为数值求解器是使用数值公差进行优化的,故严格不等式在求解时没有任何意义。因此,如果我们需要设置严格上限,则必须设置一个边界。边界的选择很重要:若其太小,则没什么作用,因为它会淹没在求解器用来定义足够接近可行约束的常规公差中;若其太大,则可能会使可行空间大大缩小。
- my_tolerance_for_strict = 1e-5;
- F = [0 <= P(1,1) <= 2-my_tolerance_for_strict];
for循环也能用于连接多个约束:
- F = [0 <= P(1,1) <= 2];
- for i = 2:n-1
- F = [F, P(i,1) <= P(2,i) - P(i,i)];
- end
在定义了变量和约束之后,就可以进行各类优化问题的求解了。接下来我们将从线性规划开始,学习一系列普通规划问题及其MATLAB实现方法。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。