YALMIP学习笔记-基础知识

7 篇文章 4 订阅

基本知识

下面的代码介绍了你可能需要的一些基础知识.YALMIP通过spdvar来定义变量,利用通过使用sdpsettings来定义约束(constraints),目标函数(objectives)以及包含设置求解器的选项(options).
通过利用optimize来求解问题,检查结果以及取得最终的解.YALMIP默认使用QUADPROG来求解问题,该优化器Matlab默认包含.

假设我们想求:
m i n      f ( x ) = x 2 + ∑ i = 1 10 ∣ x i ∣ s . t . ∑ i = 1 10 x i ≤ 10 x 1 = 0 0.5 ≤ x 2 ≤ 1.5 x 1 + x 2 ≤ x 3 + x 4 x 2 + x 3 ≤ x 4 + x 5 x 3 + x 4 ≤ x 5 + x 6 x 4 + x 5 ≤ x 6 + x 7 min\;\; f(x) = x^2 + \sum_{i=1}^{10}|x_i| \\ s.t. \sum_{i=1}^{10} x_i \leq 10\\ x_1 = 0\\ 0.5 \leq x_2 \leq 1.5\\ x_1 + x_2 \leq x_3 +x_4\\ x_2 + x_3 \leq x_4 + x_5 \\ x_3 + x_4 \leq x_5 + x_6\\ x_4 + x_5 \leq x_6 + x_7 minf(x)=x2+i=110xis.t.i=110xi10x1=00.5x21.5x1+x2x3+x4x2+x3x4+x5x3+x4x5+x6x4+x5x6+x7
则可以通过下面的方式来实现:

%定义变量
x = sdpvar(10, 1);

打印结果:

Linear matrix variable 10x1 (full, real, 10 variables)
Coeffiecient range: 1 to 1

% 定义约束
Constraints = [sum(x) <= 10, x(1)==0, 0.5 <= x(2) <= 1.5];	
for i = 1:7
	Constraints = [Constraints, x(i) + x(i+1)  <= x(i+2) + x(i+3)];
end

打印结果:

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Coefficient range|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Element-wise inequality 1x1| 1 to 10|
| #2| Equality constraint 1x1| 1 to 1|
| #3| Element-wise inequality 2x1| 0.5 to 1.5|
| #4| Element-wise inequality 1x1| 1 to 1|
| #5| Element-wise inequality 1x1| 1 to 1|
| #6| Element-wise inequality 1x1| 1 to 1|
| #7| Element-wise inequality 1x1| 1 to 1|
| #8| Element-wise inequality 1x1| 1 to 1|
| #9| Element-wise inequality 1x1| 1 to 1|
| #10| Element-wise inequality 1x1| 1 to 1|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%定义目标函数
Objective = x' * x + norm(x,1);

输出结果:

Nonlinear scalar (real, 11 variables)
Coeffiecient range: 1 to 1

%为YALMIP以及求解器设置一些选项
options = sdpsettings('verbose', 1, 'solver', 'quadprog', 'quadprog.maxiter', 100)
%求解问题
sol = optimize(Constraints, Objective, options);
%分析错误标志
if sol.problem == 0
	% 说明计算成功,此时展示得到的结果
	solution = value(x)
else
	display('错了亲');
	sol.info
	yalmiperror(sol.problem)
end

输出结果:

Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the selected value of the constraint tolerance.
solution =
0
0.5000
0.0833
0.4167
0.1667
0.3333
0.2500
0.2500
0.3333
0.1667

以上是利用YALMIP求解的一个简单例子.下面介绍基本概念.

YALMIP的标识符

YALMIP中最重要的命令是sdpvar.该命令被用来定义决策变量.
如果想要定义一个矩阵(或标量)P,该矩阵有n行m列,我们可以利用如下代码实现:

P = sdpvar(n, m);

square matrix 默认是对称的.为获得全参数化的square matrix, 可以在定义时引入第三个参数:

P = sdpvar(3,3,'full');

第三个参数可以被用来获得一些变量的预定义类型,例如Toeplitz,Hankel,diagonal, symmetric以及skew-symmetric矩阵.另外标准的Matlab语句同样适用于一个向量.

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));

YALMIP才提供查看变量的功能,所以这里的sdpvar的值普通版是看不到的.但是他会告诉你该变量的一些基本属性,以供你进行判断.
在一些情形中,可以利用多为变量来简化代码.YALMIP同样支持这种方式,可以通过cell数组或者时多为sdpvar对象来实现.

cell数组的形式如下所示:

for i = 1:5
	X{i} = sdpvar(2, 3);
end

通过在sdpvar中使用向量的维度,同样的cell数组可以被设置为如下的形式:

X = sdpvar([2 2 2 2 2],[3 3 3 3 3]);

cell数组的使用与在Matlab中的使用相同.

上面方法有个缺点,就是变量X不能直接使用,作为标准的sdpvar对象(例如加法操作等并不会在Matlab中重载).作为候选,一个完整的通用的多维sdpvar是可行的,我们可以利用下面的代码创建一个等价的对象:

X = sdpvar(2, 3, 5);

与使用cell的方式不同的是,我们可以利用标准的Matlab代码直接操作这个对象.

Y = sum(X, 3);			# 在第三个维度上对X进行求和
X(:, :, 2);					# 显示X中的第2个元素

如果前两个维度是一样的,根据标准的YALMIP的语法,可以通过下面的代码实现:

X = sdpvar(2, 2, 2, 2, 'full');

约束

我们把一系列的约束条件串联起来,从而实现约束条件的定义.

约束的意思就是上下文相关.可以通过下面的代码来实现对称矩阵以及正定约束:

n = 3;
P = sdpvar(n, n);
C = [P>= 0];

同时元素为正的正定矩阵被定义为:

P = sdpvar(n, n);
C = [P(:) >= 0];

注意到这样的定义对于非对角线上的元素约束了两次.一个好的SDP求解器将会在与处理中检测到并减少这种模式,但是同样可以利用标准的Matlab代码手动的定义独立的元素:

C = [triu(P) >= 0];

或者

C = [P(find(triu(ones(n)))) >= 0];

基于上述的规律,一个非二次型矩阵(非对称)中的正的元素,可以利用>=操作符快速得到.

P = sdpvar(n, 2*n);
C = [P >= 0];

约束可以通过简单的add以及拼接来进行定义:

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| Coefficient range|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Matrix inequality 3x3| 1 to 1|
| #2| Element-wise inequality 1x1| 1 to 2|
| #3| Equality constraint 1x1| 1 to 10|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

对于两边都有的约束,可以利用double-side来定义:

F = [0 <= P(1, 1) <= 2];

严格不等式可以是不存在的,因为各种优化器求解出来的东西都需要有一个容许的误差范围,为此我们可以定义如下形式的不等式,来等价于严格不等式(带有裕量的):

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

在定义完变量(variable)约束(constraints)之后,此时就具备了解决问题的基本条件.

优化问题的求解

线性规划

在这个例子中,我们给出两组数据,分别称之为蓝色和绿色.我们的目标就是利用线性分类器,尽可能地将这两组数据分开

blues = randn(2, 25);
greens = randn(2, 25) +	 2;

将上面地数据可视化:

plot(greens(1,:), greens(2, :), 'g*')
hold on
plot(blues(1, :) blues(2, :), 'b*')

在这里插入图片描述

线性分类器意味着我们想要寻找一组向量 a a a和标量 b b b,使得对于所有的蓝色数据点都满足 a T x + b ≥ 0 a^T x+b \geq 0 aTx+b0,且对所有的蓝色数据点都满足 a T x + b ≤ 0 a^T x+b \leq 0 aTx+b0(也就是寻找一个分离平面).通过观察数据可以发现,找到一个直线将上面的两组数据完全分开是不可能的.为此我们的目的变成了寻找一个分离平面,使得错分的数据点尽可能地少.

作为错误分类的代理, 我们引入正数 u u u v v v,并且将分类改为 a T x + b ≥ 1 − u a^Tx+b\geq 1-u aTx+b1u,且 a T x + b ≤ − ( 1 − v ) a^T x+b\leq -(1-v) aTx+b(1v).如果 u u u v v v都很小,此时我们可以得到一个较好的分离平面(margin越大).

我们定义决策变量如下:

a = sdpvar(2, 1);
b = sdpvar(1);

我们将会对每个样本点使用一个变量 u u u以及 v v v,因此我们创建了两个合适长度的向量,接下来会介绍为什么要定义成行向量的形式.

u = sdpvar(1, 25);
v = sdpvar(1, 25);

分类约束可以定义如下:

Constraints= [a '*greens+ b> = 1-u,a' *blues+ b <= -(1-v),u> = 0,v> = 0]


a T × g r e e n + b ≥ 1 − u a T × b l u e s + b ≤ − ( 1 − v ) a^T\times green + b \geq 1-u\\ a^T\times blues + b \leq -(1-v) aT×green+b1uaT×blues+b(1v)
我们想要u和v变小,在某些情况下,这意味着此时的分类性能是好的.实现这一目的的一个简单的方法就是最小化u和v的所有元素的和.然而在这种形势下条件不够充足,因此我们将绝对值的约束添加到所有a的元素中,使其所有元素都小于1.
m i n &ThickSpace;&ThickSpace; ∑ i u i + ∑ j v j s . t . &ThickSpace;   − 1 ≤ a ≤ 1 min\;\;\sum_i u_i + \sum_j v_j\\ s.t.\;\ -1 \leq a \leq 1 miniui+jvjs.t. 1a1

Objective=sum(u)+sum(v)
Constraints = [Constraints,-1 <= a <= 1];

最后,我们已经可以解决这个问题了:

optimize(Constraints, Objective)

最优的 a a a b b b通过value来获得.为更好的说明结果,可视化得到的分离平面如下所示:

x = sdpvar(2,1); 
P1 = [-5 <= X <= 5,value(a)' * x +value(b) >= 0];
P2 = [-5 <= X <= 5,value(a)' * x +value(b) <= 0];
clf
plot(P1);hold on
plot(P2);
plot(greens(1,:),greens(2,:), 'g*')
plot(blues(1,:),blues(2,:), 'b*')

遇到了错误,仅作演示使用,有发现错误原因的请告知

  • 36
    点赞
  • 249
    收藏
    觉得还不错? 一键收藏
  • 11
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 11
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值