mpt3学习记录--day2
mpt_demo_sets2
简介:Basic manipulation with convex sets(凸集的基本操作)
- YSet的创建
- YSet的存储信息
- 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');
一些例子:
- 不等式约束
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);
- 超平面约束(更一般的情况,可以是曲线)
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);
- 矩阵约束形式(高维)
% 创建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);
- 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);
- 约束相加
x = sdpvar(2,1);
% constraint1
box = ([0;0.5] <= x <= [1;2]);
% constriant2
circle = (2*norm(x-1) <=1);
S8 = YSet(x,box + circle);
- 在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
- intersect(交集)
- affine maps(仿射映射)
- comparison(比较)
- 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
- 缩放(scaling)
R = Polyhedron('lb',[1,1],'ub',[3,2]);
R1 = R*0.6; % 注意,左乘或右乘无影响,因为是数乘运算
plot(R,'color','g',R1,'color','r');
- 缩放+旋转(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');
- 平移(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');
- 使用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