We have just learnt how the logical conditions that encode the individual components of discrete hybrid automata can be turned into linear inequalities. Here we summarize them all into a single mathematical model called mixed logical dynamical (MLD).
The first two (sets of) equations resemble the state and output equations of a linear time-invariant (LTI) system, but with some additional terms corresponding to auxiliary variables (plus an offset). The auxiliary variables are not completely arbitrary – they must comply with the third (set of) equation(s).
NoteNote on the “nonsymmetry” in the format
Admittedly, our grouping of both the continuous-valued and binary-valued state variables into a single state vector, and similarly for the input and output vectors, while still keeping the auxiliary continuous-value and binary-valued variables separate, is somewhat “nonsymmetric”. Indeed, we could perhaps group the continuous-value and binary auxiliary variables into a single vector as well: \bm w(k) = \begin{bmatrix}\bm z(k) \\ \bm \delta(k) \end{bmatrix}, which will give us the alternative MLD format
\begin{aligned}
\bm x(k+1) &= \mathbf A\bm x(k) + \mathbf B_u \bm u(k) + \mathbf B_w \bm w(k) + \mathbf B_0,\\
\bm y(k) &= \mathbf C\bm x(k) + \mathbf D_u \bm u(k) + \mathbf D_w \bm w(k) + \mathbf D_0,\\
\mathbf E_w \bm w(k) &\leq \mathbf E_u \bm u(k) + \mathbf E_x \bm x(k) + \mathbf E_0,
\end{aligned}
but we decided to stick to the format supported by the main references on MLD systems.
Uniqueness of the state and output responses predicted by MLD models
An immediate question that must pop up in our minds is that of the uniqueness of the state and output responses predicted by the MLD model. For given \bm x(k) and \bm u(k), uniqueness of the values of the auxuliary variables \bm \delta(k) and \bm z(k) compliant with the inequalities, hence the next state \bm x(k+1) and the output \bm y(k), is not immediately obvious. As claimed in the seminal paper [1], MLD models corresponding to real systems do have unique solutions.
Examples
Example 1 (Switched linear system) We consider a discrete-time switched system that switches between two linear systems
x(k+1) = a_i x(k) + b_i u(k), \quad i\in\{1,2\},
where i=1 if x(k)< 0 and i=2 otherwise.
The auxiliary binary variable \delta (actually \delta(k), one for each discrete time k) is introduced to encode the switching condition:
[\delta(k) = 1] \leftrightarrow [-x(k)\leq 0].
This equivalence can be rewritten as
\begin{aligned}
-x(k) &\leq (1-\delta(k))M,\\
-x(k) &\geq \epsilon + (m-\epsilon)\delta(k),
\end{aligned}
where m and M are lower and upper bounds on x(k), respectively, and \epsilon is a small positive number. Say, if it is known that x\in [-10,10], then m=-10, M=10. We can set \epsilon = 10^{-8}.
Furthermore, just a single continuous (real) variable z is introduced to encode one of two modes of the system:
\begin{aligned}\\
[\delta(k) = 1] \rightarrow [z(k) &= a_2 x(k) + b_2 u(k)],\\
[\delta(k) = 0] \rightarrow [z(k) &= a_1 x(k) + b_1 u(k)].
\end{aligned}
This can be rewritten as
\begin{aligned}
(m_1-M_2)\delta(k) + z(k) &\leq a_1 x(k) + b_1 u(k),\\
(m_2-M_1)\delta(k) - z(k) &\leq -a_1 x(k) - b_1 u(k),\\
(m_2-M_1)(1-\delta(k)) + z(k) &\leq a_2 x(k) + b_2 u(k),\\
(m_1-M_2)(1-\delta(k)) - z(k) &\leq -a_2 x(k) - b_2 u(k),
\end{aligned}
where m_1 and M_1 are lower and upper bounds on a_1 x(k) + b_1 u(k), and m_2 and M_2 are lower and upper bounds on a_2 x(k) + b_2 u(k).
On top of these inequalities, we need to add the state update equation. It is, however, particularly simple in this case:
x(k+1) = z(k),
from which we can extract the remaining matrices parameterizing the MLD model in Eq. 1:
\begin{aligned}
A &= 0, \quad B_u = 0, \quad B_\delta = 0, \quad B_z = 1, \quad B_0 = 0,\\
C &= 1, \quad D_u = 0, \quad D_\delta = 0, \quad D_z = 0, \quad D_0 = 0.
\end{aligned}
Note, in particular, that A=0, which might be confusing at first, because it does not correspond to any of the time-discretized LTI models, but it is just a demonstration of the fact that MLD is not a state equation, it is a different type of a model.
The inequality parameterized by \mathbf E_{\delta}, \mathbf E_z, \mathbf E_x, \mathbf E_x and \mathbf E_0, for given \bm x(k) and \bm u(k), defines a set in the \delta-z plane. When the integrality of \delta is taken into consideration, the resulting set is just a singleton (a single point), which confirms the uniqueness of the state and output responses predicted by the MLD model, see Fig. 1.
Figure 1: The set of feasible values of the auxiliary variables \delta and z. The blue polyhedron relaxes the integrality condition on \delta, but if the condition is taken into consideration, the set is a single (orange) point.
Example 2 (Switched linear system with memory) Consider the first-order linear system
x(k+1) = a x(k) + b_{q(k)} u(k), \quad q(k)\in\{0,1\},
where the variable q(k) is a discrete state variable that that evolves – with some (perhaps transparent) abuse of notation – according to
q(k+1) = q(k) \lor [x(k) \leq x_{\text{lb}}],
and x_\text{lb}=-1, \; a=0.5, \; b_1=0.1, \; b_2=0.3. The initial continuous state is x(0) = 1, consequently the initial discrete state is q(k)=1. The control input u(k) is constrained to the interval [-10,10].
Obviously, once the real state x(k) falls below x_{\text{lb}}, the system switches to the other mode and stays there forever.
NoteNote on the common abuse of notation
If we wanted to avoid the aformentioned abuse of notation, we could introduce a Boolean variable Q(k) and write
Q(k+1) \leftrightarrow (Q(k) \lor [x(k) \leq x_{\text{lb}}]),
where Q(k) \leftrightarrow [q(k)=1]. Or we could write even without explicitly introducing the variable Q(k):
[q(k+1)=1] \leftrightarrow ([q(k)=1] \lor [x(k) \leq x_{\text{lb}}]),
but the somewhat abusive view of q as both a logical (Boolean) and binary (integer) variable is a common practice.
To comply with the notation of this framework, we composed the state vector of the the real and binary components as
\bm x(k) = \begin{bmatrix} x(k)\\ q(k)\end{bmatrix},
and in what follows we are going to refer to the real and binary state variables as x_\mathrm{c}(k) and x_\mathrm{b}(k), respectively.
The auxilliary binary variable (actually variables, one for each time) \delta_1(k), is introduced to encode the switching condition (the EG component of the discrete hybrid automaton):
[\delta_1(k) = 1] \leftrightarrow [x_\mathrm{c}(k) - x_\text{lb} \leq 0].
We can rewrite this as
\begin{aligned}
x_\mathrm{c}(k) - x_\text{lb} &\leq (1-\delta_1(k))M,\\
x_\mathrm{c}(k) - x_\text{lb} &\geq \epsilon + (m-\epsilon)\delta_1(k),
\end{aligned}
where m and M are lower and upper bounds on x_\mathrm{c}(k), respectively, and \epsilon is a small positive number. Say, it is assumed that x_\mathrm{c}\in [-10,10], then m=-10, and M=10. We can set \epsilon to something like \epsilon = 10^{-8}.
We can also write the mode selection (MS) component of the discrete hybrid automaton as
\text{if}\; [\delta_1(k) = 1],\; \text{then}\; z(k) = a x_\mathrm{c}(k) + b_1 u(k), \; \text{else} \; z(k) = a x_\mathrm{c}(k) + b_2 u(k).
This we can rewrite as a set of linear inequalities
\begin{aligned}
(m_2-M_1)\delta_1(k) + z(k) &\leq a x_\mathrm{c}(k) + b_2 u(k),\\
(m_1-M_2)\delta_1(k) - z(k) &\leq -a x_\mathrm{c}(k) - b_2 u(k),\\
(m_1-M_2)(1-\delta_1(k)) + z(k) &\leq a x_\mathrm{c}(k) + b_1 u(k),\\
(m_2-M_1)(1-\delta_1(k)) - z(k) &\leq -a x_\mathrm{c}(k) - b_1 u(k),
\end{aligned}
where m_1 and M_1 are lower and upper bounds on a x_\mathrm{c}(k) + b_1 u(k), and m_2 and M_2 are lower and upper bounds on a x_\mathrm{c}(k) + b_2 u(k).
We also need to encode the finite state automaton/machine (FSM) component of the discrete hybrid automaton. We write down the state update equation once again here:
[x_\mathrm{b}(k+1)=1] \leftrightarrow [x_\mathrm{b}(k)=1] \lor [\delta_1(k) = 1]
and for convenience, we rewrite it using (temporarily introduced) Boolean variables:
C \leftrightarrow A \lor B.
The equivalence can be rewritten as
C \rightarrow (A \lor B), \qquad (A\lor B) \rightarrow C.
The former can be rewritten as
\neg C \lor (A \lor B),
which (when mapping to the original binary variables) is equivalent to
(1-x_\mathrm{b}(k+1)) + x_\mathrm{b}(k) + \delta_1(k) \geq 1,
which can be finally writen as
x_\mathrm{b}(k+1) \leq x_\mathrm{b}(k) + \delta_1(k).
Note that this inequality contains the state at the next time, which is not supported by the MLD formalism. But the fix is easy. We introduce another auxilliary binary variable \delta_2(k) \coloneqq x_\mathrm{b}(k+1) and rewrite the inequality as
\delta_2(k) \leq x_\mathrm{b}(k) + \delta_1(k).
The latter implication can be rewritten as
\neg (A\lor B) \lor C,
which can be rewritten as
(\neg A\land \neg B) \lor C,
in which the distributive law gives
(\neg A \lor C) \land (\neg B \lor C),
which can finally be rewritten using the original binary variables as two inequalities
\begin{aligned}
(1-x_\mathrm{b}(k)) + x_\mathrm{b}(k+1) &\geq 1,\\
(1-\delta_1(k)) + x_\mathrm{b}(k+1) &\geq 1.
\end{aligned}
Once again, replacing the state at the next time by the auxilliary binary variable \delta_2(k) \coloneqq x_\mathrm{b}(k+1) (and subtracting the 1 from both sides), we get
\begin{aligned}
-x_\mathrm{b}(k) + \delta_2(k) &\geq 0,\\
-\delta_1(k) + \delta_2(k) &\geq 0.
\end{aligned}
We are now ready to start collecting the equations and inequalities into the MLD model. The state update equation is
\begin{bmatrix}
x_\mathrm{c}(k+1)\\
x_\mathrm{b}(k+1)
\end{bmatrix}
=
\begin{bmatrix}
0 & 0\\
0 & 1
\end{bmatrix}
\begin{bmatrix}
\delta_1(k)\\
\delta_2(k)
\end{bmatrix}
+
\begin{bmatrix}
1\\
0
\end{bmatrix}
z(k).
The inequalities are
\begin{aligned}
M\delta_1(k) &\leq -x_\mathrm{c}(k) + M + x_\text{lb},\\
(m-\epsilon)\delta_1(k) &\leq -x_\mathrm{c}(k) - \epsilon - x_\text{lb},\\
(m_2-M_1)\delta_1(k) + z(k) &\leq a x_\mathrm{c}(k) + b_2 u(k),\\
(m_1-M_2)\delta_1(k) - z(k) &\leq -a x_\mathrm{c}(k) - b_2 u(k),\\
(m_1-M_2)(1-\delta_1(k)) + z(k) &\leq a x_\mathrm{c}(k) + b_1 u(k),\\
(m_2-M_1)(1-\delta_1(k)) - z(k) &\leq -a x_\mathrm{c}(k) - b_1 u(k),\\
-\delta_1(k) + \delta_2(k) &\leq x_\mathrm{b}(k),\\
-\delta_2(k) &\leq -x_\mathrm{b}(k),\\
\delta_1(k) - \delta_2(k) &\leq 0.
\end{aligned}
Two general lessons can be learnt from the previous example. First, the procedure of finding the matrices defining the MLD model is rather tedious and error-prone. Second, the MLD model is, indeed, not a state-space model, but a different kind of a model. In particular, in the latter example we can see that the control input u does not affect the next state through the state update equation, it only appears in the inequalities.
HYSDEL language
To address the first lesson learnt from the examples – that the procedure for turning the discrete-time hybrid automaton into a MLD model is rather tedious and error-prone –, some automation of the procedure is desirable. The only tool for this that we are aware of is the HYSDEL language and parser (see the section on software for details on availability). Here we do not explain it – we refer the student to the section 16.7 in [2], which is freely downloadable –, but we show the code for Example 1, which is perhaps self-explaining.
SYSTEM switched_LTI_1 {
INTERFACE {
INPUT {
REAL u [-1.0,1.0];
}
STATE {
REAL x [-10, 10];
}
OUTPUT {
REAL y;
}
PARAMETER {
REAL A1 = -3/4;
REAL A2 = 1/5;
REAL B1 = 1;
REAL B2 = 1;
}
}
IMPLEMENTATION {
AUX {
REAL z;
BOOL delta;
}
AD {
delta = x >= 0;
}
DA {
z = {IF delta THEN A2*x + B2*u ELSE A1*x + B1*u};
}
CONTINUOUS {
x = z;
}
OUTPUT {
y = x;
}
}
}
With modern optimization modelling languages (OML), we can formulate the (optimization) problem by entering the inequalities and equations directly. It is then the responsibility of the OML to transform the problem definition into the format that the chosen MILP or MIQP solver accepts. Below we demonstrate this by means of examples.
Example 3 Here we solve the problem from Example 1 but we do not formulate the matrices from Eq. 1 explicitly. We use JuMP optimization modelling layer (for Julia language) for that. Note that JuMP also allows writing the indicators constraints directly, avoiding the need to introduce the big-M and small-m constants (but see the section on software for details for a discussion of caveats).
Figure 2: Simulation of the switched linear system from Example 1 using its MLD model without explicitly formulating the MLD matrices
Example 4 (Heating and cooling system modeled as MLD system without explicitly formulating the MLD matrices) The example is from the collection of demos for the Hybrid Toolbox for Matlab (see the section on software for details on availability). It is also described in the accompanying presentation, page 13.
We consider a room with a heater and a cooler. The temperature is sensed at two places (corresponding to workplaces of two persons). The logic of the temperature control system is this: if the person 1 feels cold, the heater is turned on; in addition, the heater can also be turned on if the person 2 feels cold, but in this case this will only happen if at the same time the person 1 does not feel hot. And similarly for the cooling: if the person one feels hot, the cooler is on; the cooler can also be on if the person 2 feels hot, but at the same time person 1 is not cold.
While the control logic is already fixed here, it is the external (ambient) temperature varying in time that serves as the input to the system. Here we fix it to a precomputed profile because our only goal at the moment is mere simulation, but it could als be optimized over when performing some verification analysis, in particular reachability analysis, but more on this later in the course.
Show the code
α1=1.0# Thermal model parameter coefficient at the location 1.α2=0.5# Thermal model parameter coefficient at the location 2.k1 =0.8# Gain of the heater/cooler at the location 1.k2 =0.4# Gain of the heater/cooler at the location 2.T1hot =30.0# Hot temperature threshold at location 1.T1cold =15.0# Cold temperature threshold at location 1.T2hot =35.0# Hot temperature threshold at location 2.T2cold =10.0# Cold temperature threshold at location 2.Q̇cmax =2.0# Air conditioning power when switched on.Q̇hmax =2.0# Heater power when switched on.Tinit = [30.0, 20.0] # Initial state (the temperatures at the two locations).tfinal =40.0# Final time.h =0.5# Sampling period.t =0:h:tfinalTa_precomp =10*cos.(t/10) .+20# Precomputed ambient temperature profile.N =length(t) # Number of time steps. Tmin, Tmax =-10.0, 50.0# Lower and upper bounds on the temperatures.ϵ =1e-6# Small constant for indicator constraints.usingJuMPusingHiGHS#using Gurobisim_problem =Model(HiGHS.Optimizer) set_silent(sim_problem)# Auxiliary binary variables (delta)@variable(sim_problem, is_1_hot[1:N], Bin) @variable(sim_problem, is_2_hot[1:N], Bin)@variable(sim_problem, is_1_cold[1:N], Bin)@variable(sim_problem, is_2_cold[1:N], Bin)@variable(sim_problem, is_heating[1:N], Bin)@variable(sim_problem, is_cooling[1:N], Bin)@variable(sim_problem, is_2_cold_and_1_not_hot[1:N], Bin)@variable(sim_problem, is_2_hot_and_1_not_cold[1:N], Bin)# Auxiliary continuous (real) variables (z)@variable(sim_problem, 0.0<= Q̇h[i=1:N] <= Q̇hmax) @variable(sim_problem, 0.0<= Q̇c[i=1:N] <= Q̇cmax) # State variables (x)@variable(sim_problem, -10.0<= T1[i=1:N+1] <=50.0) # [x(0), x(1), ..., x(N)] @variable(sim_problem, -10.0<= T2[i=1:N+1] <=50.0) # Must be literals here.# Input variables (u, ambient temperature here)@variable(sim_problem, -10.0<= Ta[i=1:N] <=50.0)for i in1:N# Event generation constraints# Is location 1 hot?@constraint(sim_problem, is_1_hot[i] --> {T1[i] >= T1hot})@constraint(sim_problem, !is_1_hot[i] --> {T1[i] <= T1hot-ϵ})# Is location 2 hot?@constraint(sim_problem, is_2_hot[i] --> {T2[i] >= T2hot})@constraint(sim_problem, !is_2_hot[i] --> {T2[i] <= T2hot-ϵ})# Is location 1 cold?@constraint(sim_problem, is_1_cold[i] --> {T1[i] <= T1cold})@constraint(sim_problem, !is_1_cold[i] --> {T1[i] >= T1cold+ϵ})# Is location 2 cold?@constraint(sim_problem, is_2_cold[i] --> {T2[i] <= T2cold})@constraint(sim_problem, !is_2_cold[i] --> {T2[i] >= T2cold+ϵ})# Is heating ON?# Is location 2 cold and location 1 not hot?@constraint(sim_problem, is_2_cold[i] - is_2_cold_and_1_not_hot[i] >=0)@constraint(sim_problem, 1- is_1_hot[i] - is_2_cold_and_1_not_hot[i] >=0)@constraint(sim_problem, -is_2_cold[i] + is_1_hot[i] + is_2_cold_and_1_not_hot[i] >=0) # Is location 1 cold or (location 2 cold and location 1 not hot)?@constraint(sim_problem, -is_heating[i] + is_1_cold[i] + is_2_cold_and_1_not_hot[i] >=0)@constraint(sim_problem, -is_1_cold[i] + is_heating[i] >=0)@constraint(sim_problem, -is_2_cold_and_1_not_hot[i] + is_heating[i] >=0)# Is cooling ON?# Is location 2 hot and location 1 not cold?@constraint(sim_problem, is_2_hot[i] - is_2_hot_and_1_not_cold[i] >=0)@constraint(sim_problem, 1- is_1_cold[i] - is_2_hot_and_1_not_cold[i] >=0)@constraint(sim_problem, -is_2_hot[i] + is_1_cold[i] + is_2_hot_and_1_not_cold[i] >=0) # Is location 1 hot or (location 2 hot and location 1 not cold)?@constraint(sim_problem, -is_cooling[i] + is_1_hot[i] + is_2_hot_and_1_not_cold[i] >=0)@constraint(sim_problem, -is_1_hot[i] + is_cooling[i] >=0)@constraint(sim_problem, -is_2_hot_and_1_not_cold[i] + is_cooling[i] >=0)# Mode-selecting constraints (here it sets just parts of the right-hand sides, not the whole z's)@constraint(sim_problem, is_heating[i] --> {Q̇h[i] == Q̇hmax})@constraint(sim_problem, !is_heating[i] --> {Q̇h[i] ==0.0})@constraint(sim_problem, is_cooling[i] --> {Q̇c[i] == Q̇cmax})@constraint(sim_problem, !is_cooling[i] --> {Q̇c[i] ==0.0})# Switched affined system dynamics@constraint(sim_problem, T1[i+1] == T1[i] + h * ( -α1*(T1[i]-Ta[i]) + k1 * (Q̇h[i]-Q̇c[i]) ) )@constraint(sim_problem, T2[i+1] == T2[i] + h * ( -α2*(T2[i]-Ta[i]) + k2 * (Q̇h[i]-Q̇c[i]) ) )# Initial condition@constraint(sim_problem, T1[1] == Tinit[1])@constraint(sim_problem, T2[1] == Tinit[2])# Precomputed ambient temperature@constraint(sim_problem, Ta[i] == Ta_precomp[i])end@objective(sim_problem, Min, 0.0)optimize!(sim_problem)
Figure 3: Simulation of the heating-cooling system using its MLD model without explicitly formulating the MLD matrices
Piecewise affine systems
Tightly related to the MLD model is another representation of a hybrid system – a PWA (piecewise affine) system (actually a model)
\begin{aligned}
\bm x(k+1) &= \mathbf A_{i(k)}\bm x(k) + \mathbf B_{i(k)} \bm u(k) + \mathbf f_{i(k)},\\
\bm y(k) &= \mathbf C_{i(k)}\bm x(k) + \mathbf D_{i(k)} \bm u(k) + \mathbf g_{i(k)},\\
& \; \mathbf H_{i(k)} \bm x(k) + \mathbf J_{i(k)} \bm u(k) \leq \mathbf K_{i(k)}.
\end{aligned}
Equivalence of MLD, DHA, and PWA (and some more)
As stated in the section 16.8 in the freely downloadable [2], the three models – MLD, DHA, and PWA – are equivalent. In fact two more kinds of models that we also covered in the course – linear complementarity systems and max-(min-)plus(-scaling) systems are also included in this equivalence statement, the key reference is [3].
A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,”Automatica, vol. 35, no. 3, pp. 407–427, Mar. 1999, doi: 10.1016/S0005-1098(98)00178-2.
W. P. M. H. Heemels, B. De Schutter, and A. Bemporad, “Equivalence of hybrid dynamical models,”Automatica, vol. 37, no. 7, pp. 1085–1091, Jul. 2001, doi: 10.1016/S0005-1098(01)00059-0.
