Return True if to_region is reachable from_region.

For details see solve_feasible.

Compute S0 \subseteq trans_set from which P2 is N-reachable.

N-reachable = reachable in horizon C{N}.
The system dynamics are C{ssys}.
The closed-loop algorithm solves for one step at a time,
which keeps the dimension of the polytopes down.

Time semantics:

- C{use_all_horizon = False}: fixed sampling period of
  discrete-valued environment variables.
  Reachability in exactly C{N} steps.

- C{use_all_horizon = True}: sampling period that varies and
  is chosen by the system, depending on how many steps are
  taken during each trajectory from C{P1} to C{P2}.
  Reachability in C{1..N} steps, with an under-approximation
  of the attractor set.

  If the system dynamics do not allow staying at the same state,
  then disconnected polytopes can arise, possibly causing issues
  in the computation. Consider decreasing the sampling period
  used for discretizing the associated continuous-time dynamical system.

@type P1: C{Polytope} or C{Region}
@type P2: C{Polytope} or C{Region}
@type ssys: L{LtiSysDyn}
@param N: horizon length
@param closed_loop: If C{True}, then take 1 step at a time.
    This keeps down polytope dimension and
    handles disturbances better.
@type closed_loop: bool

@param use_all_horizon: Used for closed loop algorithm.

    - If C{True}, then check for reachability in C{< N} steps,

    - otherwise, in exactly C{N} steps.

@type use_all_horizon: bool

@param trans_set: If specified,
    then force transitions to be in this set.
    Otherwise, P1 is used.

@return: states from which P2 is reachable
@rtype: C{Polytope} or C{Region}

Compute the components of the polytope:

   L [x(0)' u(0)' ... u(N-1)']' <= M

which stacks the following constraints:

• x(t+1) = A x(t) + B u(t) + E d(t)
• [u(k); x(k)] \in ssys.Uset for all k

If list_P is a Polytope:

• x(0) \in list_P if list_P
• x(k) \in Pk for k= 1,2, .. N-1
• x(N) \in PN

If list_P is a list of polytopes:

• x(k) \in list_P[k] for k= 0, 1 ... N

The returned polytope describes the intersection of the polytopes for all possible

Parameters:

• ssys (LtiSysDyn) - system dynamics
• N - horizon length
• disturbance_ind - list indicating which k's that disturbance should be taken into account. Default is [1,2, ... N]
• Pk (Polytope)
• list_P (list of Polytopes or Polytope)
• PN (Polytope)

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.

Parameters:

• G - The matrix to maximize with respect to
• D - Polytope describing the disturbance set
• N - Horizon length

Returns:

d_hat: Array describing the maximum possible effect from the disturbance