Proposition preserving partition module.
Partition class with following fields
domain: the domain we want to partition, type: polytope
num_prop: number of propositions
list_region: proposition preserving regions, type: list of Region
num_regions: length of the above list
adj: a matrix showing which regions are adjacent
trans[i,j] = 1 means state i is reachable from state j.
list_prop_symbol: list of symbols of propositions
orig_list_region: original proposition preserving regions
new region
is active in that region to each region in ppp
Main function that takes a domain (state_space) and a list of propositions (cont_props), and returns a proposition preserving partition of the state space.
This function takes a proposition preserving partition and generates another proposition preserving partition such that each part in the new partition is a convex polytope
This function takes a piecewise affine system pwa_sys and a proposition preserving partition ppp whose domain is a subset of the domain of pwa_sys and returns a refined proposition preserving partition where in each region a unique subsystem of pwa_sys is active.
Modified from Petter Nilsson’s code implementing merge algorithm in Nilsson et al. Temporal Logic Control of Switched Affine Systems with an Application in Fuel Balancing, ACC 2012.
Input:
Output:
Algorithms related to discretization containing both MATLAB interface and Python discretization. Input calculation is available in python only but should work also for a state space partition discretized in MATLAB.
CtsSysDyn class for specifying the continuous dynamics:
s[t+1] = A*s[t] + B*u[t] + E*d[t] + K u[t] in Uset - polytope object d[t] in Wset - polytope object
A CtsSysDyn object contains the fields A, B, E, K, Uset and Wset as defined above.
Constructor:
CtsSysDyn ([ A = [][, B = [][, E = [][, K = [][, Uset = [][, Wset = []]]]]]])
PwaSubsysDyn class for specifying a subsystem of piecewise affine continuous dynamics.
s[t+1] = A_i*s[t] + B_i*u[t] + E_i*d[t] + K_i u[t] in Uset_i - polytope object d[t] in Wset_i - polytope object for H_i*s[t]<=g_i - subdomain, type: polytope object
PwaSubsysDyn class inherits from CtsSysDyn with the additional field:
PwaSysDyn class for specifying a piecewise affine system. A PwaSysDyn object contains the fields:
For the system to be well-defined the sub_domains of its subsystems should be mutually exclusive (modulo intersections with empty interior) and cover the domain.
Compute the components of the polytope L [x(0)’ u(0)’ ... u(N-1)’]’ <= M which stacks the following constraints
If list_P is a Polytope:
If list_P is a list of polytopes
The returned polytope describes the intersection of the polytopes for all possible Input:
ssys: CtsSysDyn dynamics
N: horizon length
list_P: list of Polytopes or Polytope
Pk, PN: Polytopes
be taken into account. default is [1,2, ... N]
Refine the partition and establish transitions based on reachability analysis.
Input:
part: a PropPreservingPartition object
ssys: a CtsSysDyn or PwaSysDyn object
N: horizon length
partition.
algorithm should be used. default True.
use_mpt: if True, use MPT-based abstraction algorithm.
to stay inside starting cell. If false, safety is ensured by keeping the sequence inside the original proposition preserving cell which needs to be convex. In order to use the value false, ensure to have a convex initial partition or use prop2partconvex to postprocess the proposition preserving partition before calling discretize.
reachability analysis.
for reach- ability also in less than N steps.
transition. a value of 1 checks transitions only between neighbors, a value of 2 checks neighbors of neighbors and so on.
non-neighbors.
abs_tol: maximum volume for an “empty” polytope
verbose: level of verbosity
Output:
Load the data from MATLAB discretize implementation.
Input:
Discretize the continuous state space using MATLAB implementation.
Input:
part: a PropPreservingPartition object
ssys: a CtsSysDyn object
N: horizon length
the MATLAB implementation of discretize.
partition.
maxNumIterations: the maximum number of iterations
closed loop algorithm. For the difference between the closed loop and the open loop algorithm, see Borrelli, F. Constrained Optimal Control of Linear and Hybrid Systems, volume 290 of Lecture Notes in Control and Information Sciences. Springer. 2003.
horizon length up to probStruct.N can be used. This option is relevant only when the closed loop algorithm is used.
the reachability problem between subcells of the original partition, the cell of the original partition should be used for the safe set.
If negative, the timeout won’t be used. Note that using timeout requires MATLAB parallel computing toolbox.
in computing reachability.
verbose: level of verbosity
Generate an input file for MATLAB implementation of discretize.
Input:
Refine the partition and establish transitions based on reachability analysis.
Input:
part: a PropPreservingPartition object
ssys: a CtsSysDyn object
N: horizon length
partition.
algorithm should be used. default False.
to stay inside starting cell. If false, safety is ensured by keeping the sequence inside the original proposition preserving cell.
in reachability analysis.
for reach- ability also in less than N steps.
abs_tol: maximum volume for an “empty” polytope
Output:
Calculates the sequence u_seq such that - x(t+1) = A x(t) + B u(t) + K - x(k) in P1 for k = 0,...N - x(N) in P3 - [u(k); x(k)] in PU
and minimizes x’Rx + 2*r’x + u’Qu
Return an integer specifying in which discrete state the continuous state x0 belongs to.
Input: - x0: initial continuous state - part: PropPreservingPartition object specifying the state space partition
Output: - cellID: int specifying the discrete state in part x0 belongs to, -1 if x0 does not belong to any discrete state.
Note1: If there are overlapping partitions (i.e., x0 can belong to more than one discrete state), this just returns the first ID
Calculate an input signal sequence taking the plant from state start to state end in the partition part, such that f(x,u) = x’Rx + r’x + u’Qu + mid_weight*|xc-x(0)|_2 is minimal. xc is the chebyshev center of the final cell. If no cost parameters are given, Q = I and mid_weight=3 are used.
Input:
x0: initial continuous state
ssys: CtsSysDyn object specifying system dynamics
space partition.
start: int specifying the number of the initial state in part
end: int specifying the number of the end state in part
N: the horizon length
size (N*xdim x N*xdim). If empty, zero matrix is used.
r: cost vector for x = [x(1)’ x(2)’ .. x(N)’]’, size (N*xdim x 1)
size (N*udim x N*udim). If empty, identity matrix is used.
mid_weight: cost weight for |x(N)-xc|_2
state during execution. if False, plant is forced to stay inside the original proposition preserving cell.
been used.
make sure that the calculated input sequence is safe.
Output: - A (N x m) numpy array where row k contains u(k) for k = 0,1 ... N-1.
Note1: The same horizon length as in reachability analysis should be used in order to guarantee feasibility.
Note2: If the closed loop algorithm has been used to compute reachability the input needs to be recalculated for each time step (with decreasing horizon length). In this case only u(0) should be used as a control signal and u(1) ... u(N-1) thrown away.
Note3: The “conservative” calculation makes sure that the plants remains inside the convex hull of the starting region during execution, i.e. x(1), x(2) ... x(N-1) are in conv_hull(starting region). If the original proposition preserving partition is not convex, safety can not be guaranteed.
Calculate the array d_hat such that d_hat = max(G*DN_extreme), where DN_extreme are the vertices of the set D^N.
This is used to describe the polytope L*x <= M - G*d_hat. Calculating d_hat is equivalen to taking the intersection of the polytopes L*x <= M - G*d_i for every possible d_i in the set of extreme points to D^N.
Input:
Output: - d_hat: Array describing the maximum possible effect from disturbance
Checks if the plant remains inside P0 for time t = 1, ... N-1 and that the plant reaches P1 for time t = N. Used to test a computed input sequence. No disturbance is taken into account.
Input:
Output: - True if x(k) in P0 for k = 1, .. N-1 and x(N) in P1. False otherwise
Computes the subset x0 of `P1’ from which `P2’ is reachable in horizon `N’, with respect to system dynamics `ssys’. The closed loop algorithm solves for one step at a time, which keeps the dimension of the polytopes down.
Input:
P1: A Polytope or Region object
P2: A Polytope or Region object
ssys: A CtsSysDyn object
N: The horizon length
polytope dimension and handles disturbances better. Default: True.
allow reachability also in less than N steps.
set. If empty, P1 is used
Output: x0: A Polytope or Region object defining the set in P1 from which
P2 is reachable
A computational geometry module for polytope computations. The module can be accessed by writing
> import tulip.polytope as pc
Polytope class with following fields
representation of a polytope
representation of a polytope
array: python array in the case of a union of convex polytopes
chebXc: coordinates of chebyshev center (if calculated)
chebR: chebyshev radius (if calculated)
bbox: bounding box (if caluclated)
running reduce)
else, use A and b without modification.
Class for lists of convex polytopes
Contains the following fields:
list_poly: list of Polytope objects
list_prop: list of propositions inside region
bbox: if calculated, bounding box of region (see bounding_box)
fully dimensional
volume: if calculated, volume of region
chebXc: coordinates of maximum chebyshev center (if calculated)
chebR: maximum chebyshev radius (if calculated)
Compute the smallest hyperbox containing the polytope or region
Calculate the Chebyshev radius and center for a polytope.
If input is a region the largest Chebyshev ball is returned.
Input: poly1: A Polytope object
Output:
rc,xc: Chebyshev radius rc (float) and center xc (numpy array)
N.B., this function will return whatever it finds in attributes chebR and chbXc if not None, without (re)computing the Chebyshev ball.
Example (low dimension):
r1,x1 = cheby_ball(P, [1]) calculates the center and half the length of the longest line segment along the first coordinate axis inside polytope P
Get the dimension of a polytope or region.
Input: polyreg: Polytope or Region object
Output: dim: Dimension of input
Compute envelope of a region.
The envelope is the polytope defined by all “outer” inequalities a x < b such that {x | a x < b} intersection P = P for all polytopes P in the region. In other words we want to find all “outer” equalities of the region.
If envelope can’t be computed an empty polytope is returned
Input: polyreg: Polytope or Region abs_tol: Absolute tolerance for calculations
Output: envelope: Envelope of input
Compute the extreme points of a _bounded_ polytope
Input: - poly1: Polytope in dimension d
Output: - A (N x d) numpy array containing the N vertices of poly1
Compute the intersection between two polytopes or regions
Input: - poly1,`poly2`: Polytopes to intersect
Output: - Intersection described by a polytope
Checks if two polytopes or regions are adjacent by enlarging both slightly and checking the intersection
Input: - poly1,poly2: Polytopes or Regions to check - abs_tol: absolute tolerance - overlap: used for overlapping polytopes, functions returns
True if polytopes are neighbors OR overlap
Output: True if polytopes are adjacent, False otherwise
Check if a region is convex.
Input: reg: Region object
Output: result,envelope: result indicating if convex. if found to be
convex the envelope describing the convex polytope is returned.
Check if the description of a polytope is empty
Input: polyreg: Polytope or Region instance
Output: result: Boolean indicating whether polyreg is empty
Check if a polytope or region has inner points.
Input: - polyreg: Polytope or Region instance
Output:
otherwise.
Checks if the point p0 satisfies all the inequalities of poly1.
Input: poly1: Polytope or Region object.
Output: result: Boolean being True or False
Compute a set difference poly1 poly2 between two regions or polytopes
Input:
Output: - region: Region describing the set difference
Return N as list of bits, zero-filled to places.
E.g., given N=7, num_bin returns [1, 1, 1, 0, 0, 0, 0, 0].
Projects a polytope onto lower dimensions.
Input:
poly1: Polytope to project
dim: Dimensions on which to project
made to choose the most suitable solver.
default is 0 (be silent).
Available solvers are:
Output: - Projected polytope in lower dimension
Example: To project the polytope P onto the first three dimensions, use
>>> P_proj = projection(P, [1,2,3])
Helper function implementing “Equality set projection”. Very buggy.
Help function implementing vertex projection. Efficient in low dimensions.
Help function implementing Fourier Motzkin projection. Should work well for eliminating few dimensions.
Helper function implementing the “iterative hull” method. Works best when projecting _to_ lower dimensions.
Use quickhull to compute a convex hull.
Input: - vertices: A N x d array containing N vertices in dimension d
Output: - Polytope describing the convex hull
Removes redundant inequalities in the hyperplane representation of the polytope with the algorithm described at http://www.ifor.math.ethz.ch/~fukuda/polyfaq/node24.html by solving one LP for each facet
Warning: - nonEmptyBounded == 0 case is not tested much.
Input: poly: Polytope or Region object
Output: poly_red: Reduced Polytope or Region object
Subtract a region from a polytope
Input: - poly: polytope from which to subtract a region - reg: region which should be subtracted - abs_tol: absolute tolerance
Output: - polytope or region containing non-overlapping polytopes
Divide a region into several regions such that they are all connected.
Input: - reg1: Region object - abs_tol: Absolute tolerance
Output: List [] of connected Regions
Compute the union of polytopes or regions
Input: - polyreg1, polyreg2: polytopes or regions - check_convex: if True, look for convex unions and simplify
Output: - region of non-overlapping polytopes describing the union
Approximately compute the volume of a Polytope or Region.
A randomized algorithm is used.
Input: - polyreg: Polytope or Region
Output: - Volume of input
Functions for plotting Polytopes and Partitions. The functions can be accessed by
> from tulip.polytope.plot import *
Functions:
- get_patch
- plot
- plot_partition
- plot_trajectory
Takes a Polytope and returns a Matplotlib Patch Polytope that can be added to a plot
Example:
> # Plot Polytope objects poly1 and poly2 in the same plot
> import matplotlib.pyplot as plt
> fig = plt.figure()
> ax = fig.add_subplot(111)
> p1 = get_patch(poly1, color="blue")
> p2 = get_patch(poly2, color="yellow")
> ax.add_patch(p1)
> ax.add_patch(p2)
> ax.set_xlim(xl, xu) # Optional: set axis max/min
> ax.set_ylim(yl, yu)
> plt.show()
Plots a 2D polytope or a region using matplotlib.
Input: - poly1: Polytope or Region
Plots a 2D PropPreservingPartition object using matplotlib
Input:
ppp: A PropPreservingPartition object
Requires transitions to be stored in ppp.
each Region.
for creating custom plots.
Plots a PropPreservingPartition and the trajectory generated by x0 input sequence u_seq.
Input: