mpt_demo_sets2

简介：Basic manipulation with convex sets（凸集的基本操作）

1. YSet的创建
2. YSet的存储信息
3. Yset的性质

YSet的创建

第一步：定义集合中的变量，例如集合中变量由两个列元素，可以定义成：

x = sdpvar(2,1);

第二步：定义约束，例如区间约束可以表示为：

F1 = [-1 <= x  <= 1];

当存在多个约束时可以采用在方括号中用分号将多个约束间隔，比如下面例子中的第1、2个；也可以采用表达式相加的形式，比如下面例子中的第5个。

第三步：创建YSet

S1 = YSet(x, F1);
% 绘图
plot(S1);

上述得到的结果为：

YSet的数组操作，可以将多个YSet定义在一个数组之中，比如：

% create 1D set:
y  =sdpvar(1);
S4 = YSet(y, [-5 <= y <= 1]);
% create 2D set:
x = sdpvar(2,1);
S5 = YSet(x, [x(2).^2+x(1)+2.1*x(2) <= 1; -x(1)+0.5*x(2).^2-2.1*x(2) <=1]);
S = [S4,S5];
plot(S4, 'color', 'r', S5, 'color', 'y'); 

一些例子：

1. 不等式约束

x = sdpvar(2,1);
F2 = [1*x(1)-0.9*x(2) <= 2.5; 0.6*x(1)+1.3*x(2) >= -2.3; 0.2*x(1)+0.9*x(2) <= 4.5];
S2  = YSet(x, F2);
plot(S2);

2. 超平面约束（更一般的情况，可以是曲线）

x = sdpvar(2,1);
F3 = [x(2)>=1.5*(x(1)+5).^2-5*x(1)-20; -3.5*x(1)+2.4*x(2) <= 28];
S3 = YSet(x,F3);
plot(S3);

3. 矩阵约束形式（高维）

% 创建X为2*2的矩阵
X = sdpvar(2);
A = [1,2;-0.5,1.6];
F4 = [A'*X + X*A <= -eye(2)];
% 注意由于X此时是矩阵，表示其中的元素要用X(:)
S6 = YSet(X(:),F4);

4. create a circle centered at the point [1;2]

x = sdpvar(2,1);
F5 = [3*norm(x-[1;2])<=5];
S7 = YSet(x,F5);
plot(S7);
% 或者写成
x = sdpvar(1,2);
F5 = [3*norm(x-[1,2])<=5];
S7 = YSet(x,F5);
plot(S7);

5. 约束相加

x = sdpvar(2,1);
% constraint1
box = ([0;0.5] <= x <= [1;2]);
% constriant2
circle = (2*norm(x-1) <=1);
S8 = YSet(x,box + circle);

6. 在MPC中的应用

A = [2 -1;1 0];
B = [1;0];
nu = 1;   %number of inputs
Q = eye(2);
R = 1;
N = 2;
x0 = [2;1];
u = sdpvar(nu,N);
constraints = [];
objective = 0;
x = x0;
for k = 1:N
    x = A*x+B*u(k);
    objective = objective + norm(Q*x,1) + norm(R*u(k),1);
    constraints = [constraints, -1 <= u(k) <=1, -5 <=x <= 5];
end
S9 = YSet(u,constraints)
plot(S9)

YSet的存储信息

如果S是YSet，则其相关信息有

S.vars           % S中的变量
S.constraints    %S中的约束

以S1为例，S1.vars为：

Linear matrix variable 2x1 (full, real, 2 variables)
Values in range [-0.15643,0.98769]
Coeffiecient range: 1 to 1

S1.constraints为：

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| ID| Constraint| Coefficient range|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
| #1| Element-wise inequality 4x1| 1 to 1|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

在S中存储结构体data

S.Data.a = [1,2];
S.Data.b = {[0.1 -0.5],[0.12 -0.51],[0.11 -0.53]}
S.Data.name = 'friction';

查看YSet的性质

S.isBounded     %是否有界
S.isEmptySet    %是否空集
S.Dim           %维数（变量的个数）

mpt_demo_sets3

1. intersect（交集）
2. affine maps（仿射映射）
3. comparison（比较）
4. convex hull（凸包）

intersect（交集）

% P采用H-representation定义
P = Polyhedron([-1 -2;-0.3 0.4;0.5 -0.4;0.5 0.5],[1.7;1.4;-1.1;-0.6]);
% Q采用V-representation定义
Q = Polyhedron([-1 1;0.5 1;1 -0.4]);
% T采用V-representation定义
T = Polyhedron([-3, -1; -2, 2; 0, -2]);
% L是P和T的交集
L = P.intersect(T);
L.isEmptySet
% K是Q和的交集
K =  Q.intersect(T);
K.isEmptySet
plot(P,'color','r',Q,'color','b',T,'color','y',L,'color','g');
legend('P','Q','T','L');

affine maps（仿射映射）

scaled, rotated and shifted

1. 缩放（scaling）

R = Polyhedron('lb',[1,1],'ub',[3,2]);
R1 = R*0.6;   % 注意，左乘或右乘无影响，因为是数乘运算
plot(R,'color','g',R1,'color','r');

2. 缩放+旋转（rotation and scaling）
   左乘矩阵A

R = Polyhedron('lb',[1,1],'ub',[3,2]);
A = [1 2; 0.5 -1];
R2 = A*R;
plot(R,'color','g',R2,'color','r');

3. 平移（shifting）

R = Polyhedron('lb',[1,1],'ub',[3,2]);
R3 = R + [4;3];
R4 = R - [4;3];
plot(R,'color','g',R3,'color','r',R4,'color','b');

4. 使用affineMap
   矩阵A是的维数是n*d
   如果n<d，降维
   如果n>d，升维
   如果n=d，不变
   语法为：S = R.affineMap(A），即S是A左乘R的结果。

R = Polyhedron('lb',[1,1],'ub',[3,2]);
% S1降维
S1 = R.affineMap([1 0.5]);
% S2不变
S2 = R.affineMap([1 0.5;-0.5 2]);
% S3升维
S3 = R.affineMap([1 0.5;-0.5 2;0.4 -1.3]);

comparison（比较）

P2 == P1 Tests if two polyhedron are equal
P2 <= P1 Tests if P2 is a subset of P1
P2 >= P1 Tests if P1 is a subset of P2

convex hull（凸包）

两个维数相同的polyhedra的凸包为

P(1) = Polyhedron([-1 -1;1 0;0 1;-1 0]);
P(2) = Polyhedron([-0.5 0;2 0.5;2 -0.5]);
H = PolyUnion(P).convexHull;
plot(P(1),'color','r',P(2),'color',[0.85,0.33,0.10]');
hold on;
V1=H.V(:,1);
[V1,pos]=sort(V1);% V1是排列之后的，pos是排序后的下标
V2=H.V(pos,2);
V = [V1,V2];
V = [V;V(1,:)]
plot(V(:,1),V(:,2),'color','k','linewidth',2);

Minkowski sum（闵可夫斯基和）

P1 = Polyhedron('H',[0.5 0 1.9;1.5 -2 -0.4;-1.7 0.8 0.4;-0.5 -0.4 -1.1]);
P2 = Polyhedron('lb',[0;0],'ub',[0.5 1]);
M = P1+P2;
plot(M,'color',[1,1,1],'linewidth',2.5);
hold on;
plot(P1,'color','r',P2,'color',[0.85,0.33,0.10]');

Similarly, one can compute Pontryagin difference by using the over loaded minus

P1 = Polyhedron('lb',-[1;1],'ub',[1;1]);
P2 = Polyhedron('lb',[-0.2;-0.2],'ub',[0.2,0.2]);
N = P1-P2;
plot(P1,'color','r',N,'linewidth',2.5,P2,'color','y')

Pontryagin difference也支持数组形式的操作（相同维数）

P1 = Polyhedron([-3.5 2;-2 4;0 2;0 -2;-2 -4;-4 -2]);

%Mirrors polytope P1 round the origin and is shown in yellow
P2 = -P1;

plot(P1,'color','r',P2,'color','y');
hold on;
P3 = Polyhedron('lb',[-0.5,-0.5],'ub',[0.5,0.5]);
S = [P1,P2]-P3;
plot(P3,'color','b',S,'linewidth',2.5)

Chebyshev ball

T = Polyhedron([-1 -3;1 -0.5;-0.5 3]);
T.computeHRep
data = T.chebyCenter;
circle = antenna.Circle('Radius',data.r,'Center',data.x);

plot(T,'color','r');
hold on;
plot(circle,'color','k','linewidth',2.5);
set(gca,'XLim',[-3 3]);
set(gca,'YLim',[-3 3]);

contains

判断给定点是否在polytope中

T = Polyhedron([-1 -3;1 -0.5;-0.5 3]);
T.contains([0.1 0.2]);

methods

to see all methods on Polyhedra type

methods(Polyhedron)

for further information consult the individual help files.
e.g.
help Polyhedron/contains
doc Polyhedron/contains