Demonstration of basic usage of the geometric library
mpt_demo1
P = Polyhedron('V', randn(10,2));Plot the polytope
P.plot('color','b','alpha',0.3,'linewidth',1,'linestyle','--'); axis off;
P
Polyhedron in R^2 with representations: H-rep (irredundant) : Inequalities 4 | Equalities 0 V-rep (irredundant) : Vertices 4 | Rays 0 Functions : none
P = Polyhedron('V', randn(10,2), 'R', randn(1,2));
P.plot;
P = Polyhedron('V', randn(50,3), 'R', randn(1,3));
P.plot; axis vis3d; figure(1);
P = Polyhedron('H', [randn(30,3) ones(30,1)], 'He', [randn(1,3) 0]);
P.plot('alpha',0.3);
P = Polyhedron('lb', -rand(3,1), 'ub', rand(3,1), 'H', [randn(10,3) 1.5*ones(10,1)]);
P.plot('alpha',0.3);
P = Polyhedron(randn(10,2), ones(10,1));
P.plot('color','g'); hold on; axis square; x = 5*randn(2,1); sol = P.project(x); pplot([x sol.x]','bo'); sep = P.separate(x); v = axis; s = Polyhedron('He', sep, 'lb', [v(1);v(3)],'ub',[v(2);v(4)]); s.plot; sol = P.interiorPoint; pplot(sol.x, 'ro');
P = Polyhedron(3*randn(10,2));
Q = Polyhedron('H', [randn(10,2) ones(10,1)]);
P.plot; hold on; Q.plot;
PQ = P+Q;
PQ.plot('alpha',0.1,'color','r');
for i=1:5 P(i) = Polyhedron('V', [0 0 0;randn(1,3)]); end
plot(P,'linewidth',2); hold on; Q = Polyhedron; for i=1:5 Q = Q + P(i); end Q.plot('alpha',0.1,'color','b'); axis vis3d
x = sdpvar(2,1);
T = eye(2);
F = [x <= 0.1 ; x'*T'*T*x <= 1];
Y = YSet(x, F);
Y.plot('alpha',0.3,'color','b','grid',50); hold on; x = [-2;-2]; sol = Y.project(x); pplot([x sol.x]','bo'); sol = Y.separate(x); v = axis; s = Polyhedron('He', sol, 'lb', [v(1);v(3)],'ub',[v(2);v(4)]); s.plot('alpha',0.2); O = Y.outerApprox; O.plot('alpha',0.1,'color','b'); axis(axis*1.1);
Plotting... 35 of 50
Ps = PolyUnion;
for i=1:5 Ps.add(Polyhedron(randn(10,2)) + 5*randn(2,1)); end
for i=1:5 Ps.add(Polyhedron('H', [randn(10,2) ones(10,1)], 'He', [randn(1,2) 0]) + 5*randn(2,1)); end
for i=1:5 Z = Polyhedron; for j=1:3, Z = Z + Polyhedron([0 0;randn(1,2)]); end Ps.add(Z+5*randn(2,1)); end
Ps.plot;
T = Polyhedron('H', [randn(30,2) ones(30,1)]);
T = triangulate(T);
T.plot;
F = Polyhedron('V', 2*randn(10,2))
Polyhedron in R^2 with representations: H-rep : Unknown (call computeHRep() to compute) V-rep (redundant) : Vertices 10 | Rays 0 Functions : none
F.addFunction(Function(@(x) sin(x(1))*cos(x(2))),'func1');
F.fplot;
func = @(x) sin(x(1))*cos(x(2));
for i=1:5 F(i) = Polyhedron('V', 2*randn(10,2)) + 5*randn(2,1); end F.addFunction(Function(@(x) sin(x(1))*cos(x(2))),'wave');
Ps = PolyUnion(F);Plot the function over the set
Ps.fplot;
n = 2; m = 1; N = 5;Setup the process model and constraints
A = [1 1; 0 1];
B = [1; 0.5];
xlb = -10*ones(n,1);
xub = 10*ones(n,1);
ulb = -5*ones(m,1);
uub = 5*ones(m,1);
R = 2*eye(m);
Q = 0.2*eye(n);Formulate MPQP using Yalmip
x = sdpvar(n,N,'full');
u = sdpvar(m,N-1,'full');
cost = 0;
F = [];
for i=1:N-1 F = F + (x(:,i+1) == A*x(:,i) + B*u(:,i)); F = F + (xlb <= x(:,i) <= xub); F = F + (ulb <= u(:,i) <= uub); if i > 1 cost = cost + x(:,i)'*Q*x(:,i); end cost = cost + u(:,i)'*R*u(:,i); end
F = F + [xlb <= x(:,end) <= xub];
cost = cost + x(:,end)'*Q*x(:,end);Solve using MPQP solver Change globally the parametric QP solver
mptopt('pqpsolver','MPQP');Construct the problem
problem1 = Opt(F, cost, x(:,1), u(:));Solve
res1 = problem1.solve;
Calling mpt_mpqp_26 with default options... mpt_mpqp: 13 regionsSolve using PLCP solver Change globally the parametric QP solver
mptopt('pqpsolver','PLCP');Call problem constructor
problem2 = Opt(F, cost, x(:,1), u(:));Solve
res2 = problem2.solve;
mpt_plcp: 13 regionsPlot the partitions
plot(res1.xopt); title('PLCP'); axis tight;
plot(res2.xopt); title('MPQP'); axis tight;
cost = 0;
F = [];
for i=1:N-1 F = F + (x(:,i+1) == A*x(:,i) + B*u(:,i)); F = F + (xlb <= x(:,i) <= xub); F = F + (ulb <= u(:,i) <= uub); if i > 1 cost = cost + norm(Q*x(:,i),1); end cost = cost + norm(R*u(:,i),1); end
F = F + [xlb <= x(:,end) <= xub];
cost = cost + norm(Q*x(:,end),1);Solve using MPLP solver
mptopt('plpsolver','MPLP');
problem3 = Opt(F, cost, x(:,1), u(:));
res3 = problem3.solve;
Calling mpt_mplp_26 with default options... mpt_mplp: 28 regionsSolve using PLCP solver
mptopt('plpsolver','PLCP');
problem4 = Opt(F, cost, x(:,1), u(:));
res4 = problem4.solve;
regions: 21, unexplored: 5 mpt_plcp: 28 regionsPlot the partitions
plot(res3.xopt); title('MPLP'); axis tight;
plot(res4.xopt); title('PLCP'); axis tight;
◀ | mpt_demo2 | mpt_demo_unions2 | ▶ |
© 2010-2013 Colin Neil Jones: EPF Lausanne, colin.jones@epfl.ch