1

Pulse-Based Control Using Koopman Operator
Under Parametric Uncertainty

arXiv:1708.00232v1 [] 1 Aug 2017

Aivar Sootla and Damien Ernst

Abstract—In applications, such as biomedicine and systems/synthetic biology, technical limitations in actuation complicate implementation of time-varying control signals. In order
to alleviate some of these limitations, it may be desirable to
derive simple control policies, such as step functions with fixed
magnitude and length (or temporal pulses). In this technical note,
we further develop a recently proposed pulse-based solution to
the convergence problem, i.e., minimizing the convergence time
to the target exponentially stable equilibrium, for monotone
systems. In particular, we extend this solution to monotone
systems with parametric uncertainty. Our solutions also provide
worst-case estimates on convergence times. Furthermore, we
indicate how our tools can be used for a class of non-monotone
systems, and more importantly how these tools can be extended to
other control problems. We illustrate our approach on switching
under parametric uncertainty and regulation around a saddle
point problems in a genetic toggle switch system.
Index Terms—Koopman Operator, Optimal Control, Isostables, Monotone Systems, Bistable Systems, Genetic Toggle Switch

I. I NTRODUCTION
Recent developments in molecular biology, neuroscience,
bioengineering among others, opened control theory to new
fields of applications. If the implementation of optimal control
methods in biological and biomedical applications is technically feasible, it may still be costly, labor intensive and/or
invasive. This drawback in part stems from actuating using
a time-varying feedback control signal, which is typically
computed using an optimal control framework such as model
predictive control, for example. Hence it may be desirable to
use (piecewise) constant control signals in order to simplify
the actuation. We motivate our formulation through concrete
limitations and constraints by considering synthetic biology as
a motivating application, but our results can be applied to a
broader class of problems.
Synthetic biology aims to engineer and control biological
functions in living cells and is an emerging field of science
with numerous applications [1]. Recently, control theoretic
regulation of protein levels in microbes was shown to be
possible (cf. [2], [3]). One of the simplest, but nevertheless important, control problems is a switching problem. For example,
it may be desirable to switch off production of a certain protein
in a cell by an external stimulus (by adding a chemical, heat
shock, exposure to light etc). Switching between two modes
Aivar Sootla is with the Department of Engineering Science,
University of Oxford, Parks Road, Oxford, OX1 3PJ, UK
aivar.sootla@eng.ox.ac.uk. Damien Ernst is with the Montefiore
Institute, University of Liège, B-4000, Belgium dernst@ulg.ac.be.
This work was performed while A. Sootla held a postdoctoral position at
University of Liège funded by F.R.S.–FNRS. Currently, A. Sootla is supported
by the EPSRC Grant EP/M002454/1.

(production of a protein is “on” and “off”) can mathematically
be expressed as switching between stable equilibria of a model.
This problem can be solved using standard optimal control
tools, however, in general, actuation has some limitations.
If, for instance, we use the addition of a chemical as an
actuation technique, then lowering the chemical concentration
in a culture can be done through dilution, but it is more labor
intensive than increasing the concentration. It may, therefore,
be desirable to derive (piecewise-)constant (or at least easy-toimplement) control policies, which can also solve the problem.
In [4], [5], it has been proposed to solve the switching
problem using temporal pulses of a fixed length τ and a fixed
magnitude µ for monotone systems (cf. [6]). Since monotone
systems are met often in biological applications [7], restricting
the class of system to monotone is not unreasonable. In [8],
the authors extended this framework by solving a convergence
problem using temporal pulses, i.e. driving the system from an
arbitrary initial point to the target exponentially stable equilibrium. This allowed a closed-loop formulation of the pulsebased control, thus providing robustness through feedback.
Furthermore, the problem was reduced to a static optimization
problem using the Koopman operator framework (cf. [9]).
In this technical note, we further develop the results of [4]
and [8]. In particular, we extend the framework to control
under parametric uncertainty using the Koopman operator and
monotone system theory. We also use our framework as a
building block to solve other control problems, for example,
we apply our method to the problem of regulation around a
saddle point.
The rest of the note is organized as follows. In Section II,
we cover the essential background on monotone systems and
Koopman operator. In Section III, we discuss how to solve the
switching problem in a closed loop setting using the Koopman
operator, while generalizing and extending the problem formulation. In Section III-B, we consider the convergence problem
under parametric uncertainty and illustrate the framework on
numerical examples. The proofs of our technical results are
found in Appendix A.
II. P RELIMINARIES
We consider a system in the following form
ẋ = f (x, p, u),

x(0) = x0 ,

(1)

where f : D × P × U → Rn , u : R≥0 → U, where
D ⊆ Rn , P ⊆ Rk , and u belongs to the space U∞ of
Lebesgue measurable functions with values from U ⊂ R
(i.e., the system is single-input). We define the flow map
φ : R × D × P × U∞ → Rn , where φ(t, x0 , p, u) is a solution

2

to the system (1) with an initial condition x0 , parameter value
p and a control signal u. If u = 0, then we call the system (1)
unforced. Throughout the note, we assume that f ∈ C 2 and we
denote the Jacobian matrix of f (x, p, 0) as J(x, p). For every
p, we assume that the eigenvectors of J(x∗ (p), p) are linearly
independent (i.e., J(x∗ (p), p) are diagonalizable), where x∗ (p)
is a stable equilibrium of ẋ = f (x, p, 0). We denote the
eigenvalues of J(x∗ (p), p) by λi (p), we enumerate them
counting their algebraic multiplicities, i.e., i = 1, . . . , n. We
also assume that the eigenvalues λj (p) are ordered according
to their real parts, i.e., <(λi (p)) ≥ <(λj (p)), for i ≤ j. If
λ1 (p) is such that <(λ1 (p)) > <(λ2 (p)), then λ1 (p) is called
the dominant eigenvalue. We will drop the dependence on p
and write x∗ , J(x∗ ), λi when the parameter value is clear
from the context.
Koopman Operator. Spectral properties of nonlinear dynamical systems ẋ = f (x) can be described through the
Koopman operator. We limit our study of the Koopman
operator to systems with f ∈ C 2 and admitting an exponentially stable equilibrium x∗ , i.e., J(x∗ ) is a Hurwitz matrix
(its eigenvalues are such that <(λj ) < 0 for all j). We
∗
∗
assume that
 B(x ) is the basin of attraction of x , that is,
B(x∗ ) = x ∈ Rn limt→∞ φ(t, x) = x∗ .
The Koopman semigroup of operators associated with ẋ =
f (x) acts on functions g : Rn → C (also called observables)
and is defined as
U t g(x) = g(φ(t, x)),

(2)

where φ(t, x) is a flow of the system, and x ∈ B(x∗ ). The
Koopman semigroup is linear [9], and we can define the Koopman eigenfunctions as nontrivial functions s : B(x∗ ) → C
satisfying U t s(x) = s(φ(t, x)) = s(x) eλt for all t > 0,
where eλt belong to the point spectrum of U t and λ’s are
called Koopman eigenvalues. Under the assumptions above, it
can be shown that the eigenvalues λj of the Jacobian matrix
J(x∗ ) are also Koopman eigenvalues (cf. [10]). Furthermore,
there exist n eigenfunctions sj ∈ C 1 (B(x∗ )) associated with
eigenvalues λj [10] and
(f (x))T ∇sj (x) = λj sj (x).

(3)

The eigenfunction s1 (x) can be computed at a given point
using Laplace averages [11]. Furthermore, the eigenfunctions
sj (x) can be estimated using linear algebraic methods [10].
An important concept related to the eigenfunction s1 , which
is associated with the Koopman eigenvalue λ1 , is called
isostables.
Definition 1: Suppose s1 (x) is a C 1 (B(x∗ )) eigenfunction corresponding to the dominant eigenvalue λ1 with
0 > <(λ1 ) > <(λj ) for j ≥ 2. The isostables
∂Bα for α ≥ 0 are boundaries of the sublevel sets s1
∗
(i.e., Bα = {x ∈ B(x
n
o )||s1 (x)| ≤ α}), or ∂Bα =
x ∈ B(x∗ ) |s1 (x)| = α . We will also use the following
n
o
notation ∂+ Bα = x ∈ B(x∗ ) s1 (x) = α and ∂− Bα =
n
o
x ∈ B(x∗ ) s1 (x) = −α for positive α.
It can be also shown that the eigenfunction s1 (x) captures
the dominant behavior of the system, and Bα → B as α → ∞.

Furthermore, trajectories with initial conditions on the same
isostable ∂Bα1 reach other isostables ∂Bα2 after a time
 
α1
1
ln
, if α2 < α1
(4)
T =
|<(λ1 )|
α2
For an eigenvalue λ1 with multiplicity µ1 , there exist µ1
eigenfunctions, and hence the isostables are not unique. However, we will consider specific eigenfunctions and isostables
in the case of monotone systems. For more information about
the isostables see [11]. We will also use the following concept
introduced in [8].
S
Definition 2: Let r : Rn ×R≥0 ×R≥0 → C {∞} be called
a pulse control function and defined as:
r(x, µ, τ ) = s1 (φ(τ, x, µ)),
where s1 is a dominant eigenfunction on B(x∗ ), and we use
the convention that r(x, µ, τ ) = ∞, if φ(τ, x, µ) 6∈ B(x∗ ).
Monotone Systems and their Spectral Properties. We define
a partial order  as follows: x  y if and only if x−y ∈ Rn≥0 .
We will also write x  y if x  y and x 6= y, and x  y
if x − y ∈ Rn>0 . We define a partial order on the space of
signals u ∈ U∞ : u  v if u(t) − v(t) ∈ R≥0 for all t ≥ 0. Let
[x, y] denote an order-interval, that is [x, y] = {z ∈ Rn |x 
z  y}. We call a set A order-convex if for all x, y in A
the interval [x, y] is a subset of A. We note that, in general,
partial orders can be induced by other cones in Rn (cf. [12]).
Definition 3: Consider a function W : Rn → R ∪
{−∞, +∞}. The effective domain of W is defined as
dom(W ) = {x ∈ Rn ||W (x)| < ∞}. We call W increasing if
W (x) ≤ W (y) for all x  y on dom(W ).
Definition 4: The system (1) is called monotone if
φ(t, x, p, u)  φ(t, y, q, v) for all t ≥ 0, and for all x  y,
p  q, and u  v.
A certificate for monotonicity involves conditions on the
vector field of the system and is usually called Kamke-Müller
conditions (cf. [6]). We finally mention spectral properties of
monotone systems ( [13]).
Proposition 1: Consider the system ẋ = f (x) with f ∈ C 2 ,
which admits an exponentially stable equilibrium x∗ with a
basin of attraction B(x∗ ). Assume that J(x∗ ) is diagonalizable
with eigenvalues λj such that <(λ1 ) > <(λj ) for all j ≥ 2.
If the system is monotone, then λ1 is real and negative, v1 ,
which is a right eigenvector of J(x∗ ) corresponding to λ1 , and
s1 , which is a Koopman eigenfunction corresponding to λ1 ,
can be chosen such that v1  0, ∇s1 (x)  0 for all x ∈ B(x∗ )
and v1T ∇s1 (x∗ ) = 1.
Without loss of generality, we assume that for monotone
systems s1 (x) is such that ∇s1 (x)  0 for all x ∈ B(x∗ ).
III. P ULSE -BASED C ONTROL U SING KOOPMAN
O PERATOR
A. Convergence to an Isostable
The main motivating problem in [8] is a convergence
problem from an arbitrary initial point x0 to an ε-ball around
the target equilibrium x• (e.g., in the Euclidean metric). This
problem can be solved using standard optimal control methods,
however, the optimal control policy can have a complicated

3

dependence on time (in the case of the open-loop control)
or state (in the case of the closed-loop control). Therefore
the authors proceeded in a different direction and formulated
another problem as follows.
Problem 1: Converging to an isostable. Consider the system
ẋ = f (x, u) at the initial state x0 and control signals
(
1 0≤t≤τ,
u(t) = µh(t, τ ), where h(t, τ ) =
(5)
0 t>τ.
Compute µ, τ such that the flow φ(t, x0 , µh(·, τ )) reaches
Bε (x• ) with ε  1 in minimum time Tconv subject to an
energy budget kukL1 = µ · τ ≤ Emax .
The authors argued that if Tconv , τ , µ obtained from
Problem 1 are such that Tconv  τ , or µ is sufficiently small,
then we get a good approximation of the convergence time
to an ε-ball around x• . In order to simplify the problem the
following assumptions are made:
A1. Let f (x, u) be C 2 in (x, u) on D × U.
A2. Let the system ẋ = f (x, 0) have an exponentially stable
equilibrium x• with a basin of attraction B(x• ) ⊆ D.
A3. Let the system be monotone on D × U.
A4. The dominant eigenfunction s1 (x) of the system ẋ =
f (x, 0) defined on the basin of attraction B(x• ) is such
that ∇s1 (x)  0 for all x ∈ B(x• ).
A5. f (x, µ1 )  f (x, µ2 ) for all x ∈ D and µ1 > µ2 ≥ 0.
Assumption A1 guarantees existence and uniqueness of
solutions of ẋ = f (x, u) and existence of C 1 eigenfunctions
around an exponentially stable equilibrium. We ensure that the
target equilibrium is exponentially stable in Assumption A2.
Note that in the case of exponentially stable monotone systems
with f ∈ C 2 (Assumptions A1-A3), we have that ∇s1 (x)  0
and f (x, µ1 )  f (x, µ2 ) for all x ∈ B(x• ), µ1 ≥ µ2 (see [6]).
Therefore Assumption A4 holds if, for example, the Jacobian
of f (x, 0) is irreducible for all x, while Assumption A5 additionally requires that f (x, µ1 ) 6= f (x, µ2 ) for all x ∈ B(x• ),
µ1 > µ2 . Therefore by making Assumptions A4 and A5 we do
not add significant restrictions in comparison to the previous
assumptions, however, we add a certain degree of regularity.
We clarify practical implications of these assumptions in what
follows. We note that monotonicity (Assumption A3) is the
most restrictive and crucial assumption of this framework. We
also remark that the space of control signals (i.e., temporal
pulses (5)) is rich enough to solve some control problems
under Assumptions A1 – A5. Due to space limitations we
refer the reader to [8] for detail.
At the first glance, the problem of convergence to Bε (x• )
(i.e., Problem 1) does not seem to be easier than the problem
of convergence to the ε-ball around x• (e.g., {x ∈ Rn |kx −
x• k2 ≤ ε}). As shown in [8], however, Problem 1 can be
solved using the following static optimization program under
Assumptions A1 – A5:
γ∗ =

min

µ≥0,τ ≥0

ln |r(x, µ, τ )| + |λ1 |τ

subject to: r(x, µ, τ ) ≤ −ε, µ · τ ≤ Emax ,

(6)
(7)

where Tconv = |λ11 | (γ ∗ − ln(ε)). It can be shown that the
function ln |r(x, µ, τ )|+|λ1 |τ is decreasing in µ and τ . Hence

the minimum is attained only while one of the constraints
is activated. Therefore, we spend all the energy budget Emax
and/or reach the isostable Bε . More specifically, if the curves
{x ∈ Rn |s1 (x) = −ε} and µ = Emax /τ intersect, then the
optimal solution lies on this intersection. If these curves do
not intersect then we need to solve:


Emax
(8)
, τ + |λ1 |τ.
γ ∗ = min ln r x,
τ ≥0
τ
The computation of the function r can be performed using
the methods to compute s1 [8], which may not be numerically
cheap for large-scale and/or stiff systems. However, solving
an optimal control problem using dynamic programming or
maximum principle generally requires (at least) a comparable
computational effort. Furthermore, we need to solve a dynamic
optimization problem instead of the static one as in (6).
In order to apply this framework to wider class of problems,
we briefly comment on the main properties of the solution
in [8]. Consider the following auxiliary control problem:
V (z, µ, β) =

inf

τ,

τ,u∈U∞ ([0,µ])

(9)

subject to ẋ = f (x, u), x(0) = z,
x(τ ) ∈ Cβ = {x̃ ∈ Rn |s1 (x̃) = β},
where β ∈ R. Under Assumptions A1 – A5, it can be shown
that the optimal solution to (9) has the following form:
(
µ if s1 (x) < β
u(x) =
,
(10)
0 if s1 (x) ≥ β
The function r(x, µ, τ ) is related to the function V (x, µ, β),
however, this relation is not bijective. The relevant properties
of r are listed below.
Proposition 2 ( [8]): Let the system ẋ = f (x, u) satisfy
Assumptions A1–A5, let s1 (x) < β and (x, µ, τ ) ∈ dom(r).
(i) r is a C 1 function on dom(r). Furthermore, ∂τ r(x, µ, τ ) >
λ1 r(x, µ, τ ), ∂µ r(x, µ, τ ) > 0 and ∇x r(x, µ, τ )  0.
(ii) If r(x, µ, τ ) ≤ 0, then ∂τ r(x, µ, τ ) > 0.
(iii) If f (x, ν)  0, then ∂τ r(x, µ, τ ) > 0 for all finite τ > 0,
and µ > ν.
Additionally, if a solution to (9) exists then V (x, µ, β) = τ
if and only if r(x, µ, τ ) = β, provided that (ii) or (iii)
holds. If Assumptions A4 and A5 do not hold, then it
can be shown that all the inequalities are not strict, that is
∇x r  0, ∂τ r(x, µ, τ ) ≥ λ1 r(x, µ, τ ), ∂µ r(x, µ, τ ) ≥ 0 in
(i), ∂τ r(x, µ, τ ) ≥ 0 in (ii) and (iii). Furthermore, in this case
if a solution to (9) exists then V (x, µ, β) = τ if and only if
τ = argminξ≥0 r(x, µ, ξ) = β.
B. Convergence Under Parametric Uncertainty
We consider the system (1) satisfying assumptions A1–
A5 for every p ∈ [p1 , p2 ]. Furthermore, we assume that
the system (1) is monotone on D × [p1 , p2 ] × R≥0 . This
is an additional assumption with respect to Assumption A3,
since we require monotonicity of the flow with respect to
parameter changes. Let s1 (x, p) denote the eigenfunction of
the system (1) for a fixed parameter p defined on the basin

4

V (z, µ, β, p) =

inf

τ,

τ,u∈U∞ ([0,µ])

(11)

subject to (1), x(0) = z, x(τ ) ∈ Aβ ,
where Aβ = {x̃ ∈ Rn |g(x̃) = β} and g ∈ C 1 and ∇g  0
on dom(g). We have the following result, the proof of which
can be found in Appendix A.
Proposition 3: Let the system ẋ = f (x, u) satisfy Assumptions A1–A5 for every parameter p, and let ∇g(x)  0.
If g(z) < β (respectively, g(z) ≥ β), then u0 (·) = µ
(respectively, u0 (·) = 0) is an optimal solution of (11).
Now provided that we fix the function g(z) for all p ∈ [p1 , p2 ]
we have the following result with the proof in Appendix A.
Proposition 4: Let the system (1) satisfy Assumptions A1–
A5 for parameter values p1 , p2 and be monotone on D ×
[p1 , p2 ]×R≥0 . Let p1  p2 , y  z, ν ≤ µ, max(g(z), g(y)) ≤
min(α, β) then:
(i) α ≤ β implies that V (z, µ, α, p2 ) ≤ V (y, ν, β, p1 ), and the
result holds with strict inequalities;
(ii) V (z, µ, α, p2 ) ≥ V (y, ν, β, p1 ) implies that α ≥ β.
The main application of this result is control of systems
under parametric uncertainty. That is, we assume that the
bounds on parameter variations are known (for example,
p ∈ [p1 , p2 ]), however, the exact values of parameters (and
hence the eigenfunctions) are unknown. In this case, we have:
V (z, µ, β, p1 ) ≥ V (z, µ, β, p) ≥ V (z, µ, β, p2 ),
for all p ∈ [p1 , p2 ], if g(z) < β. Hence, given an initial condition, target set, and the maximum allowed control, we can
estimate the optimal value function. Note that in this case the
function g is not necessarily an eigenfunction for (1) with the
parameter p. But letting g(z) being equal to an eigenfunction
s1 (z, p1 ) (respectively, s1 (z, p2 )) we can estimate the upper
bound (respectively, the lower bound) on the value function.
In order to do so, we simply compute the function r as in the
non-parametric case.
Estimation of the convergence time to the equilibrium is
slightly more delicate. We use an extended definition of the
function r:
r(x, µ, τ, p) = s1 (φ(τ, x, p, µ), p).
The convergence time T can be computed in the same manner,
provided that the parameter value p is given:


|r(x, µ, τ, p)|
1
ln
.
T (x, µ, τ, p, ε) =
|λ1 (p)|
ε
In the non-parametric case, it was shown that ∂µ r(x, µ, τ ) > 0
and ∂τ r(x, µ, τ ) > 0 (if r(x, µ, τ ) < 0) under Assumptions A1 – A5. Therefore, if we increase µ or τ we can
predict the changes in r(x, µ, τ ) and, as a consequence, in
T (x, µ, τ, ε). In the parameter dependent case, estimating the
dependence of the function r(x, µ, τ, p) on p is not straightforward. For example, ∂r(x, µ, τ, p)/∂p does not necessarily
have the same sign pattern for fixed x, µ, and τ , while varying p. Furthermore, λ1 (p) now also depends on parameters.
However, since we are dealing with monotone systems, it is

Time to reach

Time to reach:

20

Pulse magnitude

of attraction B(x• (p)) of the exponentially stable equilibrium
x• (p). We first formulate an auxiliary problem:

15

10

With parameter
values:

5

2

4

6

Pulse length

8

Fig. 1. The level sets of T = |λ 1(q)| ln
1
qmin , qint , and qmax , where ε = 10−14 .



10

|r(x∗ (q),µ,τ,q)|
ε



for q equal to

still possible to bound the convergence time provided that the
functions r(·, ·, ·, p1 ) and r(·, ·, ·, p2 ) are computed.
Proposition 5: Let the system (1) satisfy Assumptions A1–
A5 for parameter values p ∈ [p1 , p2 ] and be monotone on
D × [p1 , p2 ] × R≥0 . Let the pair (µ, τ ) belong to the set
n
Sσ (x, p1 , p2 , ε) = (ν, ξ) ∈ R2 r(x, ν, ξ, p1 ) ≥ −εe|λ1 (p1 )|σ ,
o
r(x, ν, ξ, p2 ) ≤ −εe|λ1 (p2 )|σ ,
for a positive ε. Then the flow φ(t, x, p, µh(·, τ )) at time σ +τ
belongs to the set Atarget = {x ∈ Rn |z1  x  z2 , z1 ∈
∂− Bε (x• (p1 )), z2 ∈ ∂− Bε (x• (p2 ))}.
The proof in is Appendix A. We note that Sσ (x, p, p, ε) =
{(ν, ξ) ∈ R2 T (x, ν, ξ, p, ε) = σ}, provided that r is negative
on this level set. In this case, for p ∈ [p1 , p2 ] we have that

Sσ (x, p1 , p2 , ε) = (µ, τ ) ∈ R2 ν1 ≤ µ ≤ ν2 , ξ1 ≤ τ ≤ ξ2
(ν1 , ξ1 ) ∈ Sσ (x, p1 , p1 , ε), (ν2 , ξ2 ) ∈ Sσ (x, p2 , p2 , ε) ,
which implies that we can compute the level sets
T (x, µ, τ, p1 , ε) and T (x, µ, τ, p2 , ε) in order to estimate
Sσ (x, p1 , p2 , ε). Note however, that these levels sets can potentially intersect for large values of µ, which is clearly seen for
large values of ε. Therefore, the level sets T (x, µ, τ, p, ε) = σ
do not lie between the level sets T (x, µ, τ, p1 , ε) = σ and
T (x, µ, τ, p2 , ε) = σ provided that p ∈ [p1 , p2 ]. Furthermore,
only for small ε, the target set Atarget provides a good
estimate of the interval [x• (p1 ), x• (p2 )], which contains x• (p).
The approximation of the convergence time can hence be
conservative for large intervals [p1 , p2 ] and for large µ.
We also remark that these results can potentially be extended
to a class non-monotone systems ẋ = f (x, u) by using
bounding monotone systems approach (cf. [14]). If there are
monotone systems ẋ = r(x, u), ẋ = g(x, u) on D × R such
that r(x, u)  f (x, u)  g(x, u) for all (x, u) ∈ D × R,
then the convergence time for the system ẋ = f (x, u) can
potentially be bounded by the convergence time of the systems
ẋ = r(x, u), ẋ = g(x, u).
C. A Case Study
As a case study we consider a genetic toggle switch system
which is met in many applications (cf. [15]). We consider the

5

The states xi represent the concentration of proteins, whose
mutual repression is modeled via a rational function. The
parameters p11 and p21 model the basal synthesis rate of
each protein. The parameters p14 and p24 are degradation rate
constants, p12 , p22 describe the strength of mutual repression,
while p13 , p23 are called Hill coefficients. We consider the
following parameter values:


q1 q2 4 1
p=
,
q3 q4 3 2
with the set of admissible parameters is equal
 to [qmin , qmax ],
1.8 950 1.2 1050 and qmax =
where qmin =
2.2 1050 0.7 950 , and
 4the order q is induced by
−1
Q = diag 1 1 −1
R≥0 . We also consider qint =

2 1000 1 1000 , where qint ∈ [qmin , qmax ]. We can
2
readily check that the model is monotone on R2≥0 ×R4≥0 ×R
 2≥0
2
with respect to diag 1 −1 R≥0 × Q × diag 1 −1 R≥0 .
For every q ∈ [qmin , qmax ] the system is bistable with two
exponentially stable equilibria x∗ and x• . In this setting, the
equilibrium x∗ has the state x2 “switched on” (x2 is much
larger than x1 ), while x• has the state x1 “switched on” (x1
is much larger than x2 ). The rest of assumptions also hold,
which can be verified by direct calculation.
First we illustrate the application of Proposition 5. We
assume that u2 = 0 and we compute u1 in order to switch
from x∗ to x• .
Problem 2. Switching Under Parametric Uncertainty. Consider a system (12) with q ∈ [qmin , qmax ]. Compute
an open-loop policy u1 = µh(·, τ ) such that the flow
φ(t, x∗ (q), q, µh(·, τ ), 0) converges to x• (q) in approximately
T + τ time units for any q ∈ [qmin , qmax ].
In Figure 1 we depict the level sets of the function
T (µ, τ, q, ε) = |λ11(q)| ln(r(x∗ (q), µ, τ, q)/ε) for parameter
values qmin , qmax , and qint . For ε = 10−14 , we can bound the
location of the level sets of T (µ, τ, qint , ε) by the level sets
of T (µ, τ, qmax , ε) and T (µ, τ, qmin , ε) thus solving Problem
2 under parametric uncertainty. The function T appears to
change monotonically in q, however, for ε = 10−2 the level
sets of T actually intersect, which rules out monotonicity of T .
Similar observations were made on different sets of parameter
values. However, despite this conservatism a set of admissible
µ, τ solving Problem 2 can still be found.
Next, we consider a slightly more complicated problem and
illustrate the application of Proposition 4.
Problem 3. Event-Based Regulation Around a Saddle Point.
Consider a system (12). Let x† (q) denote the saddle point, and
let x† (q) ∈ [z, y]. Keep the flow φ(t, x, q, u1 , u2 ) in [z, y] for
any q ∈ [qmin , qmax ] for all t ≥ T for a large enough T > 0.
First we note that qmin q q q qmax implies that
B(x• (qmin )) ⊆ B(x• (q)) ⊆ B(x• (qmax )), and B(x∗ (qmin )) ⊇
B(x∗ (q)) ⊇ B(x∗ (qmax )) according to results in [16], which
is illustrated in the left panel of Figure 2. Given the geometry

Basins of attraction

10

10

9

9

8

8
7

7

6

6

5

5
4

Controlled Trajectory

State

State

following model with pij ≥ 0 for all i = 1, 2, j = 1, . . . , 4:
p12
ẋ1 = p11 +
− p14 x1 + u1 ,
1 + xp213
(12)
p22
ẋ2 = p21 +
p23 − p24 x2 + u2 .
1 + x1

4

5

6

7

State

8

9

10

4
4

5

6

7

State

8

9

10

Fig. 2. In the left panel, boundaries of basins of attraction of the toggle switch
for different parameter values. In the right panel a trajectory controlled by an
event-based procedure in order to keep the flow inside the box [4, 10]×[4, 10].

of the systems, we derived an ad-hoc event-based solution as
follows.
Proposed Event-Based Control Scheme: Let δ  1.
If x = φ(t, x, p, u1 , 0) 6∈ [z, y] compute u2 so that
φ(t, x, p, 0, u2 ) reaches ∂− B1/δ (x• (qmin )) in minimum-time.
If x = φ(t, x, p, 0, u2 ) 6∈ [z, y] compute u1 so that
φ(t, x, p, u1 , 0) reaches ∂+ B1/δ (x∗ (qmax )) in minimum time.

T

T
Let y = 10, 4 , z = 4, 10 . First, we consider the
treatment of the event φ(t, x, p, u1 , 0) 6∈ [z, y]. We pick a finite
number of points x̂i on the boundary of [z, y], and compute
the magnitudes for the control signals u2 = µh(·, ξupper ) with
ξupper = 10 steering the flow from x̂i to ∂− B1/δ (x• (qmin )).
Note that we will reach ∂− B1/δ (x• (q)) in time τ < 10
due to Proposition 4. For a small enough δ we also have
that B1/δ (x• (qmin )) ⊆ B1/δ (x• (qmax )) around the saddle
points. Therefore we can also estimate the lower bound on τ
by computing ξlower such that r(x̂i , µ, ξlower , qmin ) = −1/δ.
However, these estimates are not accurate for a large δ since
isostables do not change monotonically under parameter variation (see [16]). When we detect that x = φ(t, x, p, u1 , 0) 6∈
[z, y] we find the closest to x points x̂i , and steer the flow from
x̂i with u2 = µh(·, ξupper ). The event x = φ(t, x, p, 0, u2 ) 6∈
[z, y] is treated in a similar way. In this case, we use u1 and
aim to reach the isostable ∂+ B1/δ (x∗ (qmax )) in 10 time units.
A resulting trajectory for qint is depicted in Figure 2. We
note the flow cannot be stabilized by a single control signal
in the set [z, y], since for all constant control signal values
only saddle points are contained in [z, y]. Therefore, limiting
the flow to a set around an unstable equilibrium is the only
available solution.
IV. C ONCLUSION
In this technical note, we studied pulse-based control approach to the convergence problem. We extended this approach
to the convergence problem for monotone systems under
parametric uncertainty, where we showed how to estimate
the convergence time to the target equilibrium. We illustrated
our results on numerical examples, and we showed how our
algorithm can be applied to different control problems.
R EFERENCES
[1] P. Purnick and R. Weiss, “The second wave of synthetic biology: from
modules to systems,” Nat. Rev. Mol. Cell Biol., vol. 10, no. 6, pp. 410–
422, jun 2009.

6

[2] A. Milias-Argeitis, S. Summers, J. Stewart-Ornstein, I. Zuleta, D. Pincus,
H. El-Samad, M. Khammash, and J. Lygeros, “In silico feedback for in
vivo regulation of a gene expression circuit,” Nat biotechnol, vol. 29,
no. 12, pp. 1114–1116, 2011.
[3] J. Uhlendorf, A. Miermont, T. Delaveau, G. Charvin, F. Fages, S. Bottani, G. Batt, and P. Hersen, “Long-term model predictive control of gene
expression at the population and single-cell levels,” Proc. Nat. Academy
Sciences, vol. 109, no. 35, pp. 14 271–14 276, 2012.
[4] A. Sootla, D. Oyarzún, D. Angeli, and G.-B. Stan, “Shaping pulses to control bistable systems analysis, computation and
counterexamples,” Automatica, vol. 63, pp. 254–264, Jan. 2016,
http://arxiv.org/abs/1409.6150v3.
[5] A. Sootla and A. Mauroy, “Shaping pulses to control monotone bistable
systems using Koopman operator,” in In Proc Symposium Nonlinear
Control Systems, Aug 2016, pp. 710–715.
[6] D. Angeli and E. Sontag, “Monotone control systems,” IEEE Trans
Autom Control, vol. 48, no. 10, pp. 1684–1698, 2003.
[7] E. Sontag, “Monotone and near-monotone biochemical networks,” Syst
Synthetic Biol, vol. 1, no. 2, pp. 59–87, 2007.
[8] A. Sootla, A. Mauroy, and D. Ernst, “An optimal control formulation
of pulse-based control using Koopman operator,” July 2016. [Online].
Available: http://arxiv.org/abs/1707.08462
[9] I. Mezić, “Spectral properties of dynamical systems, model reduction
and decompositions,” Nonlinear Dynam, vol. 41, no. 1-3, pp. 309–325,
2005.
[10] A. Mauroy and I. Mezić, “Global stability analysis using the eigenfunctions of the koopman operator,” IEEE Trans Autom Control, vol. 61,
no. 11, pp. 3356–3369, 2016.
[11] A. Mauroy, I. Mezić, and J. Moehlis, “Isostables, isochrons, and koopman spectrum for the action–angle representation of stable fixed point
dynamics,” Physica D, vol. 261, pp. 19–30, 2013.
[12] H. Smith, Monotone dynamical systems: an introduction to the theory
of competitive and cooperative systems. Am Math Soc., 2008, vol. 41.
[13] A. Sootla and A. Mauroy, “Geometric properties and computation of
isostables and basins of attraction of monotone systems,” To appear in
IEEE Trans Autom Control, 2017, https://arxiv.org/pdf/1705.02853.
[14] J.-L. Gouzé, A. Rapaport, and M. Z. Hadj-Sadok, “Interval observers
for uncertain biological systems,” Ecological modelling, vol. 133, no. 1,
pp. 45–56, 2000.
[15] T. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetic
toggle switch in escherichia coli,” Nature, vol. 403, pp. 339–342, 2000.
[16] A. Sootla and A. Mauroy, “Properties of isostables and basins of
attraction of monotone systems,” in Proc Am Control Conf, 2016, pp.
7365–7370.

ν, z  y, p2  p1 , which entails that g(φ(t, z, µ, p2 )) ≥
g(φ(t, y, ν, p1 )). Hence g(φ(ξ, z, µ, p2 )) ≥ g(φ(ξ, y, ν, p1 )) =
β, which implies that the flow φ(t, z, µ, p2 ) has entered the set
Cβ− at some time t ≤ ξ, so that V (z, µ, β, p2 ) ≤ V (y, ν, β, p1 ).
Additionally since Cα− ⊇ Cβ− we have that
V (z, µ, β, p2 ) ≥ V (z, µ, α, p2 ), for β ≥ α,

(14)

since it takes more time to reach Cβ− with the same optimal
control signal. Now using (13) we get V (y, ν, β, p1 ) ≥
V (z, µ, β, p2 ) ≥ V (z, µ, α, p2 ) and complete the proof.
If α < β, then using the same arguments we can obtain the following chain of inequalities V (y, ν, β, p1 ) ≥
V (z, µ, β, p2 ) > V (z, µ, α, p2 ) for µ ≥ ν, z  y, p2  p1
(ii) We prove this part by contradiction. Let V (z, µ, α, p2 ) ≥
V (y, ν, β, p1 ) and α < β for some µ ≥ ν, z  y . According
to the point (i), we have that V (z, µ, α, p2 ) < V (y, ν, β, p1 ),
which contradicts our assumption.
Proof of Proposition 5: Let (µ, τ ) belong to Sσ (x, p1 , p2 , ε),
and let s1 (φ(τ, x, p1 , µ), p1 )
≥
−εe−λ1 (p1 )T and
−λ1 (p2 )σ
s1 (φ(τ, x, p2 , µ), p2 ) ≤ −εe
. Due to monotonicity
for all p ∈ [p1 , p2 ] we have
φ(τ + σ, x, p1 , µh(·, τ ))  φ(τ + σ, x, p, µh(·, τ )) 
φ(τ + σ, x, p2 , µh(·, τ )) .
Now, since s1 (·, p) is increasing and s1 (φ(t, x, p, 0), p) =
s1 (x, p)eλ1 (p)t for any p ∈ [p1 , p2 ], we have
s1 (φ(τ + σ, x, p, µh(·, τ )), p1 ) ≥
s1 (φ(τ + σ, x, p1 , µh(·, τ )), p1 ) =
s1 (φ(τ, x, p1 , µ), p1 )eλ1 (p1 )σ ≥ −ε.

(15)

Similarly, we can obtain, the following chain of inequalities:
s1 (φ(τ + σ, x, p, µh(·, τ )), p2 ) ≤

A PPENDIX A
P ROOFS

s1 (φ(τ + σ, x, p2 , µh(·, τ )), p2 ) =

Proof of Proposition 3: This proof is along the lines of the
proof of Theorem 8 in [8], therefore we only provide a sketch.
Let s1 (z) < β, while the case s1 (z) ≥ β is shown similarly.
Due to the premise of Proposition 3, s1 (z1 ) ≥ s1 (z2 ) for
all z1  z2 . Using Assumptions A3 – A5, it can be shown
that φ(t, x, µ)  φ(t, x, u) for all t > 0 and any u(t) such that
µ ≥ u(t) for all t > 0, therefore g(φ(τ, x, µ)) ≥ g(φ(τ, x, u)).
Therefore, by using the control signal u0 (t) = µ for all t
we reach the level set Aβ at least as fast as we any other
control signal u(t) ≤ µ for all t ≥ 0. Since we consider only
u ∈ U∞ ([0, µ]) the control signal u0 (·) = µ is an optimal
solution to the considered problem.
Proof of Proposition 4: (i) Let Cα− = {x|g(x) ≥ α}, Cβ− =
{x|g(x) ≥ β}. Due to this definition we have that Cα− ⊇ Cβ−
if and only if α ≤ β. By the premise g(z) < α, g(y) < β,
hence we can view the problem as reaching the sets Cα− , Cβ− .
We proceed by showing that
V (y, ν, β, p1 ) ≥ V (z, µ, β, p2 ),

(13)

for all z  y, µ ≥ ν, p2  p1 . Due to monotonicity we
have that φ(t, z, µ, p2 )  φ(t, y, ν, p1 ) for all t ≥ 0, µ ≥

s1 (φ(τ, x, p2 , µ), p2 )eλ1 (p2 )σ ≤ −ε.

(16)

By noticing that the inequalities (15), (16) are the definition
of the set Atarget , we complete the proof.

