mpt_demo1

Purpose

Demonstration of basic usage of the geometric library

Syntax

mpt_demo1

Description

Basic usage of the new interface to geometric library

Example(s)

Example 1

Create polytope:
P = Polyhedron('V', randn(10,2)); 
 Plot the polytope 
 P.plot('color','b','alpha',0.3,'linewidth',1,'linestyle','--'); axis off; 

../../../fig/mpt/demos/mpt_demo1_img_1.png

Double description created and stored automatically
 P 
Polyhedron in R^2 with representations:
    H-rep (irredundant) : Inequalities   4 | Equalities   0
    V-rep (irredundant) : Vertices   4 | Rays   0
Functions : none

Example 2

Create Polyhedron
P = Polyhedron('V', randn(10,2), 'R', randn(1,2)); 

P.plot;

../../../fig/mpt/demos/mpt_demo1_img_2.png

Create another Polyhedron
 P = Polyhedron('V', randn(50,3), 'R', randn(1,3)); 

 P.plot;  axis vis3d; figure(1); 

../../../fig/mpt/demos/mpt_demo1_img_3.png

Example 3

Lower-dimensional polyhedra
 P = Polyhedron('H', [randn(30,3) ones(30,1)], 'He', [randn(1,3) 0]); 

 P.plot('alpha',0.3); 

../../../fig/mpt/demos/mpt_demo1_img_4.png

Example 4

Boxes
 P = Polyhedron('lb', -rand(3,1), 'ub', rand(3,1), 'H', [randn(10,3) 1.5*ones(10,1)]); 

 P.plot('alpha',0.3); 

../../../fig/mpt/demos/mpt_demo1_img_5.png

Example 5

Polyhedron queries Project onto polyhedron:
 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'); 
            

../../../fig/mpt/demos/mpt_demo1_img_6.png

Addition
 P = Polyhedron(3*randn(10,2)); 

 Q = Polyhedron('H', [randn(10,2) ones(10,1)]); 

 P.plot; hold on; Q.plot; 

../../../fig/mpt/demos/mpt_demo1_img_7.png

 PQ = P+Q; 

 PQ.plot('alpha',0.1,'color','r'); 

../../../fig/mpt/demos/mpt_demo1_img_8.png

Lower-dimensional addition
                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
            

../../../fig/mpt/demos/mpt_demo1_img_9.png

Example 6

Operations with convex sets
 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

../../../fig/mpt/demos/mpt_demo1_img_10.png

Example 7

Create set of polytopes
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; 

../../../fig/mpt/demos/mpt_demo1_img_11.png

Create complex
 T = Polyhedron('H', [randn(30,2) ones(30,1)]); 

 T = triangulate(T); 

 T.plot; 

../../../fig/mpt/demos/mpt_demo1_img_12.png

Example 8

Polyhedral functions
 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; 

../../../fig/mpt/demos/mpt_demo1_img_13.png

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

../../../fig/mpt/demos/mpt_demo1_img_14.png

Example 9

Parametric solutions Formulate MPC problem with the following dimensions
 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 regions
Solve 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 regions
Plot the partitions
 plot(res1.xopt);  title('PLCP'); axis tight; 

../../../fig/mpt/demos/mpt_demo1_img_15.png

 plot(res2.xopt);  title('MPQP'); axis tight; 

../../../fig/mpt/demos/mpt_demo1_img_16.png

Formulate MPLP problem in YALMIP
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 regions
Solve using PLCP solver
 mptopt('plpsolver','PLCP'); 

 problem4 = Opt(F, cost, x(:,1), u(:)); 

 res4 = problem4.solve; 
regions:   21, unexplored: 5 
mpt_plcp: 28 regions
Plot the partitions
 plot(res3.xopt);  title('MPLP'); axis tight; 

../../../fig/mpt/demos/mpt_demo1_img_17.png

 plot(res4.xopt);  title('PLCP'); axis tight; 

../../../fig/mpt/demos/mpt_demo1_img_18.png

See Also

mpt_demo_sets1, mpt_demo_functions1, mpt_demo_unions1


© 2010-2013 Colin Neil Jones: EPF Lausanne, colin.jones@epfl.ch