An Optimal Control Formulation of Pulse-Based Control Using Koopman Operator
Aivar Sootla Alexandre Mauroy Damien Ernst
math.OC cs.SY

Direct Observation of Orbital Hybridisation in a Cuprate Superconductor 
    J. Chang
    December 30, 2023
========================================================================




	In many applications, and in systems/synthetic biology in particular, it is desirable to compute control policies that force the trajectory of a bistable system from one equilibrium (the initial point) to another equilibrium (the target point), or in other words to solve the switching problem. It was recently shown that, for monotone bistable systems, this problem admits easy-to-implement open-loop solutions in terms of temporal pulses (i.e., step functions of fixed length and fixed magnitude). In this paper, we develop this idea further and formulate a problem of convergence to an equilibrium from an arbitrary initial point. We show that this problem can be solved using a static optimization problem in the case of monotone systems.	Changing the initial point to an arbitrary state allows building closed-loop, event-based or open-loop policies for the switching/convergence problems. In our derivations, we exploit the Koopman operator, which offers a linear infinite-dimensional representation of an autonomous nonlinear system. One of the main advantages of using the Koopman operator is the powerful computational tools developed for this framework. Besides the presence of numerical solutions, the switching/convergence problem can also serve as a building block for solving more complicated control problems and can potentially be applied to non-monotone systems. We illustrate this argument on the problem of synchronizing cardiac cells by defibrillation. Potentially, our approach can be extended to problems with different parametrizations of control signals since the only fundamental limitation is the finite time application of the control signal.


Monotone Systems, Koopman Operator, Isostables, Generalized Repressilator, Genetic Toggle Switch




§ INTRODUCTION


Synthetic biology is an active field of research with applications in metabolic engineering, bioremediation and energy sector (<cit.>). One of the main goals of synthetic biology is to engineer biological functions in living cells (<cit.>) and control theory naturally plays an essential role toward that end. Control theoretic regulation of protein levels in microbes was achieved by <cit.>, <cit.>, and <cit.>. However, these works result in time-varying feedback control signals, which are affected by limitations due to physical constraints not only in sensing but also in actuation. Actuation limitations are observed with chemical induction, for instance. While the concentration of a chemical can be easily increased by adding this chemical to a culture, it is much more complicated to decrease its concentration (for instance, diluting is a labor-intensive process that is not suited to be performed repeatedly). Therefore, one of the challenges in synthetic biology is to derive control policies that not only achieve the desired objectives but are also simple enough to deal with the actuation limitations. In this context, an example of such control policies is based on temporal pulses, on which the present paper focuses.

One of the basic but nevertheless important control problems is the problem of convergence to an equilibrium. In synthetic biology, it can, for instance, be motivated by the genetic toggle switch circuit, which is a major building block in applications. The genetic toggle switch by <cit.>, for example, consists of two interacting genes. Their design ensured that the concentration of proteins expressed by one gene is always much higher than the concentration of proteins expressed by the other gene: one gene is switched “on” while the other is switched “off”. In this case, the system is bistable and the control objective of the problem is to drive the state from one equilibrium (e.g. one gene switched “on”) to the other (e.g. the other gene switched “on”) in minimum time and given a fixed energy budget (e.g. to avoid cell death). In this paper, we propose to solve this convergence problem using a specific set of signals, namely  
temporal pulses u(t) with fixed length τ and magnitude μ: 

    u(t) = μ h(t,τ)      
    h(t,τ) = 
     1     0 ≤ t ≤τ ,
    
    0     t >τ .

The main challenge in solving this convergence problem is the absence of well-developed theory and computational methods. This can be explained by the fact that classical optimal control methods rely on the semigroup property (of the value function and/or the dual variable), while our parametrization of control signals violates it. We, therefore, restrict our analysis to the class of monotone systems, properties of which allow simpler solutions. We note, however, that monotone systems play an important role in systems/synthetic biology and hence there exist many applications with this restriction. We solve the convergence problem by computing only one time-independent function, which we call the pulse control function and denote by r. This function links all the tunable parameters of the convergence problem so that finding a tradeoff between the convergence time and the energy budget becomes straightforward.

In our solution, we use the Koopman operator (cf. <cit.>), which offers a linear infinite dimensional description of a nonlinear system and allows a spectral approach to the nonlinear system. In particular, the definition of the pulse control function involves the dominant eigenfunction of the Koopman operator of the unforced system (i.e., when the control signal is equal to zero). This allows to use computational methods developed in the Koopman operator framework to solve our problem. In particular, we show that the function r can be computed with Laplace averages (<cit.>).

Our theoretical results generally do not apply to non-monotone systems. Let alone a solution using pulses may not exist. However, we can still use the function r in some cases to derive control policies for non-monotone systems. Furthermore, we can use our framework to solve more complicated control problems. We illustrate this idea on the problem of  synchronization of cardiac cells modeled by non-monotone FitzHugh-Nagumo systems (<cit.>).

Related Work. <cit.> proposed to solve the problem of switching between exponentially stable equilibria in monotone systems using temporal pulses. The authors derived a computational procedure to estimate the set of control signals allowing the switch. The control strategy proposed by <cit.> is open loop so that the control signal cannot be adjusted during the experiment. <cit.> considered a similar setting, but using the Koopman operator framework to estimate the set of all switching pulses and provide estimates for convergence time. In this paper, we present an optimization program for solving the convergence problem, which is a generalization of the switching problem considered by <cit.> and <cit.>.
	
<cit.> considered the convergence/escape problem cast in the Koopman framework formalism. <cit.> did not restrict the class of systems but assumed full controllability, i.e. control signals affect all the states in an affine manner. <cit.> proposed to synchronize the cardiac cells by using geometric tools from Koopman operator framework and the techniques based on the Pontryagin's maximum principle. The authors did not parametrize the control signal, which led to complicated time-dependent optimal control signals. We will show that our scheme still achieves synchronization, but in the longer time than the optimal solution by <cit.>.

To summarize, in this paper, we consider a convergence problem similarly to <cit.>, the target set is chosen similarly to <cit.>, while the class of systems and the control signals are restricted as in the works by <cit.> and <cit.>. We opted for restricting the class of systems in order to show optimality of our approach in this specific case, which indicates that systems' properties can be used to derive easy-to-implement and interpretable optimal solutions. 

Outline of the Paper. The rest of the paper is organized as follows. In Section <ref>, we cover basic definitions and properties of monotone systems and we introduce the Koopman operator framework. We formulate and discuss our convergence problem in Section <ref>, while also presenting the main result, which we prove in Section <ref>. In Section <ref>, we apply the theoretical results to the switching problem (open and closed loop control) and to the synchronization of cardiac cells.



§ PRELIMINARIES



Consider a system of the form

    ẋ = f(x,u),   x(0) = x_0,
 
with f: ×→^n, u:→, and where ⊂^n, ⊂ are open and u belongs to the space _∞ of Lebesgue measurable functions with values from . We assume that f(x, u) is twice continuously differentiable (C^2) in (x, u) on ×. 
The flow map ϕ: ××_∞→^n induced by the system is such that ϕ(t, x_0, u) is a solution of (<ref>) with an initial condition x_0 and a control signal u. We denote the differential of a function g(x, y):^n×^m→^k with respect to x as ∂_x g(x, y).  Let J(x) denote the Jacobian matrix of f(x,0) (i.e., J(x) = ∂_x f(x,0)/∂ x).  
For every stable fixed point x^∗ of the system, we assume that the eigenvectors of J(x^∗) are linearly independent (i.e., J(x^∗) is diagonalizable). The eigenvalues of J(x^∗) are denoted by λ_i with i = 1,…, n and are ordered by their real part, that is (λ_i) ≥(λ_j) for all i, j.  
We also denote the positive orthant by ^n = {x∈^n | x_i > 0,  i = 1,…, n}, and the nonnegative orthant by ^n = {x∈^n | x_i ≥ 0,  i = 1,…, n}.



 §.§ Koopman operator


Autonomous nonlinear systems can be studied in the framework of the Koopman operator. A semigroup of Koopman operators acts on functions g:^n→ (also called observables) and is defined by

    U^t g(x) = g ∘ϕ(t,x,0),      t ≥ 0

where ∘ is the composition of functions. Provided that the vector field and observables are C^1, one can define the infinitesimal generator of the operator as L g(x) = (f(x, 0))^T ∇ g(x) on a compact set. The semigroup is linear (cf. <cit.>) and can be studied through its spectral properties. In this paper, we will limit our use of the Koopman operator to unforced systems (<ref>) on a basin of attraction of an exponentially stable equilibrium x^∗ (that is, the eigenvalues λ_j of J(x^∗) are such that (λ_j) < 0 for all j). The basin of attraction is defined by (x^∗) = {x∈^n | lim_t→∞ϕ(t,x,0) = x^∗}. In this case, the eigenvalues λ_j of the Jacobian matrix J(x^∗) are also the so-called Koopman eigenvalues, which are associated with the Koopman eigenfunctions s_j: → satisfying

    U^t s_j(x) = s_j(ϕ(t, x,0)) = s_j(x)   e^λ_j t,    x∈,

or equivalently

    ∇ s_j(x)^T f(x, 0) = λ_j s_j(x).

If the vector field f(·, 0) is a C^2 function and the Jacobian matrix J(x^∗) is diagonalizable, then the Koopman eigenfunctions s_j belong to C^1 (<cit.>).

We refer to an eigenvalue λ_1 satisfying (λ_1) > (λ_j) for all λ_1 λ_j as the dominant eigenvalue. We assume that such an eigenvalue exists (it is the case for monotone systems) and we call the associated eigenfunction s_1 the dominant eigenfunction. If the dominant eigenvalue is such that (λ_1)<0, then the dominant eigenfunction s_1 can be computed through the Laplace average

    g_λ^∗(x) = lim_t→∞1/T∫_0^T (g∘ϕ(t, x,0)) e^-λ t d t.

For all g∈ C^1 that satisfy g(x^∗)=0 and (∇ g(x^∗))^T v_1 ≠ 0, where v_1 is the right eigenvector of J(x^∗) corresponding to λ_1, the Laplace average g_λ_1^∗ is equal to s_1(x) up to a multiplication with a scalar. If the algebraic and geometric multiplicities of λ_1 are equal to μ_1, then there are μ_1 independent eigenfunctions associated with λ_1 and they can be computed by choosing linearly independent right eigenvectors of J(x^∗) corresponding to λ_1. The eigenfunctions s_j(x) with j≥ 2 are generally harder to compute using Laplace averages, but can be obtained through other methods such as linear algebraic methods by <cit.>. The eigenfunctions can also be estimated from data by using the so-called dynamic mode decomposition (DMD) method (cf. <cit.>). 

The Koopman eigenfunctions capture important geometric properties of the system. In particular, the dominant Koopman eigenfunction s_1 is related to the notion of isostables. 

 Suppose that s_1 ∈ C^1. An isostable ∂_α associated with the value α > 0 is the boundary of the set
_α =  {x∈^n  |  |s_1(x)| ≤α}, that is

    ∂_α =  {x ∈^n  |  |s_1(x)| = α}.



A more rigorous definition based on the flow is due to <cit.>. Isostables are useful from a control perspective since they capture the dominant (or asymptotic) behavior of the unforced system. Indeed, it can be shown that trajectories with initial conditions on the same isostable ∂_α_1 converge synchronously toward the equilibrium, and reach other isostables ∂_α_2, with α_2<α_1, after a time

    = 1/|(λ_1)|ln(α_1/α_2) .

In the case λ_1 ∈, for example, it can be shown that the trajectories starting from ∂_α share the same asymptotic evolution

    ϕ(t,x,0) → x^∗ + v_1  α  e^λ_1 t ,    t→∞ .




 §.§ Monotone Systems

We consider systems that preserve a partial order induced by a nonnegative orthant ^n, but our results can be extended to other cones  in ^n. 
We define a partial order ≽ as follows: x≽ y if and only if x - y ∈^n (we write x⋡y if the relation  x ≽ y does not hold). We will also write x≻ y if x≽ y and x y, and x≫ y if x- y ∈^n.  Similarly, a partial order can be defined on the space of signals u∈_∞: u≽ v if u(t) - v(t) ∈^n for all t≥ 0. 
We also introduce concepts that are important for our subsequent discussion. Let [x, y] denote an interval in the order ≼, that is [x, y] = { z ∈^n | x≼ z ≼ y }. 	
For a function W:^n →∪{-∞, +∞}, we refer to  the set (W) = {x∈^n | |W(x)| < ∞} as its effective domain. A function W:^n →∪{-∞, +∞} is called increasing if W(x) ≥ W(y) for all x≽ y on (W). Control systems in the form (<ref>) whose flows preserve a partial order relation ≽ are called monotone systems. 

The system ẋ = f(x,u) is called monotone if ϕ(t,x, u)≼ϕ(t,y, v) for all t≥ 0, and for all x≼ y, u≼ v. 

 The unforced system ẋ = f(x, 0) is strongly monotone 
if it is monotone and x≺ y implies that ϕ(t,x, 0)≪ϕ(t,y, 0) for all t> 0.

A certificate for monotonicity is a condition on the vector field, for which we refer the reader to <cit.>. We finally consider the spectral properties of unforced monotone systems that are summarized in the following result. 
 
Consider that the system ẋ = f(x) with f∈ C^2() has an exponentially stable equilibrium x^∗ and let λ_j be the eigenvalues of J(x^∗) such that (λ_i)  ≥(λ_j) for all i ≤ j.

(i)  If the system is monotone with respect to  on a set ⊆(x^∗), then λ_1 is real and negative, the right eigenvector v_1 of J(x^∗) can be chosen such that v_1 ≻ 0, while the eigenfunction s_1 can be chosen such that s_1(x) ≥ s_1(y) for all x, y∈ satisfying x≽ y. 

(ii) Furthermore, if the system is strongly monotone with respect to  on a set ⊆(x^∗) then λ_1 is simple, real and negative, λ_1 > (λ_j) for all j ≥ 2, v_1 and s_1 can be chosen such that v_1≫ 0 and s_1(x)> s_1(y) for all x, y∈ satisfying x≻ y; 

The result is from <cit.>. Without loss of generality, we will assume that a dominant eigenfunction s_1 is increasing even if λ_1 is not simple. 



§ CONVERGENCE TO AN ISOSTABLE PROBLEM
 


 §.§ Problem Formulation and Discussion

In order to formulate the basic problem we want to address, consider the following assumptions:

	
  A1. The vector field f(x,u) in (<ref>) is twice continuously differentiable in (x, u) on ×. 
	
  A2. The unforced system (<ref>) has an exponentially stable equilibrium x^∗ in  with a diagonalizable J(x^∗). 
	
  A3. The system is monotone with respect to ^n× and forward-invariant on × in the sense that for all x∈, u∈_∞, the flow ϕ(t,x,u(·)) belongs to . 
	
  A4. The eigenfunction s_1(x) is such that 
∇ s_1(x) ≫ 0 for all x∈(s_1). 
	
  A5. f(x, μ_1) ≻ f(x, μ_2) for all x∈ and μ_1> μ_2≥ 0.
	
  A6. The space of control signals is limited to temporal pulses u(t) = μ h(t,τ), where h is defined in (<ref>).
	
Assumption A1 guarantees existence and uniqueness of solutions, while Assumption A2 introduces a reference point x^∗. These assumptions are perhaps more restrictive than the ones usually met in control theory. That is, f(x,u) is usually assumed to be Lipschitz continuous in x for every fixed u, and the equilibria are asymptotically stable. Our assumptions are guided by our consequent use of the Koopman operator. Assumptions A1 and A2 guarantee the existence of continuously-differentiable eigenfunctions on the basin of attraction (x^∗) of x^∗. Monotonicity is crucial, but forward-invariance on × is a rather technical assumption on which our computational methods do not rely. Assumption A4 is well-posed since s_1∈ C^1 due to Assumption A1. If Assumptions A1 – A3 hold, then we have ∇ s_1(x) ≻ 0 and f(x,μ_1) ≽ f(x, μ_2) for μ_1 > μ_2 ≥0, hence Assumptions A4 and A5 serve as technical assumptions that guarantee uniqueness of solutions and a certain degree of regularity. We will comment throughout the paper on the case when Assumptions A4 – A5 do not hold. Assumption A6 is guided by many applications, where there is a need to parametrize in advance the control signal. In this paper, we choose the easiest parametrization, although the only fundamental limitation is to set u(t) to zero after some time τ. It is, therefore, possible to generalize our approach to more complicated control signals. We proceed by formulating a basic but fundamental problem.

Problem 1. Converging to an isostable. Consider the system ẋ = f(x,u) satisfying assumptions A1–A6 and the initial state x_0. Compute a control signal u(t) = μ h(t,τ) such that the flow ϕ(t,x_0, u(·)) reaches the set _ε(x^∗) for some small ε>0 in minimum time units _ conv subject to the energy budget u_≤_ max.

Our formulation based on the isostables is guided by our use of Koopman operator framework for computational purposes. However, there are other benefits in this formulation. One can view _ε(x^∗) as a ball in the (contracting) pseudometric d_K(x,y) = |s_1(x) - s_1(y)| on a basin of attraction (x^∗)  (d_K(x, y) is a pseudometric, since d_K(x, y) can be equal to zero for some x y). By reformulating the standard convergence problem using a pseudometric defined through Koopman eigenfunctions, we take into account the dynamical properties of the unforced system. Furthermore, if ε is close to zero then the solution of Problem 1 can be used to solve a convergence-type problem. For example, the problem of switching between equilibria, which is considered in <cit.>, falls into this category. 

The main challenge in solving this problem is the parametrization of the control signal. Most of the control methods (such as dynamic programming, Pontryagin's maximum principle) are not tailored to deal with time parametrized control signals since they rely on the semigroup property of the value function or the dual variable. Hence it is not entirely clear how to systematically approach this problem through these methods.



 §.§ A Solution using a Static Optimization Program

The computational solution to our problem will be established by computing first a static function, which we call the pulse control function and define below. 

  Let the function r:××→⋃{∞} such that 
	
    r(x, μ,τ) = s_1(ϕ(τ, x, μ)),

	where s_1 is a dominant eigenfunction on the basin of attraction of x^∗, be called the pulse control function. By convention r(x, μ, τ) = ∞, if ϕ(τ, x, μ) ∉(x^∗)


If s_1 is real-valued and increasing on (s_1) = (x^∗), then we assume it is extended to ^n so that s_1:^n→∪{-∞, +∞} is increasing on ^n. We note that the switching function proposed in <cit.> corresponds to r(x^∗, ·,·). The pulse control function r can be used in a context broader than the switching function by <cit.>, which we demonstrate in this paper. In particular, we will solve Problem 1 using the following result, which we will prove in what follows.



	Consider the system (<ref>), Problem 1 under Assumptions A1–A6 and the optimization program:

    γ^∗= min_μ≥ 0, τ≥ 0     1/|λ_1|  ln(|r(x, μ,τ)|) + τ, 
    subject to:       r(x,μ,τ)≤ -ε, 
       μ·τ≤_ max.
	
	If s_1(x_0)≤-ε, an optimal solution to (<ref>) is an optimal solution to Problem 1, if the former is feasible. Furthermore, the objective is nonincreasing in μ and τ and an optimal solution to (<ref>), if it exists, is achieved at the boundary of the admissible set to the constraint (<ref>) and/or the constraint (<ref>).


If Assumptions A5 and A6 do not hold, then we can possibly have multiple minima including the points which do not activate the constraints. However, we can still compute a minimizing solution using the same program. Intuitively, the objective function is the convergence time (in fact, _ conv=γ^∗-1/|λ_1| ln(ε)), the constraint on r(x, μ, τ) ≤ - ε ensures that we stop when reaching _ε and the constraint μ·τ≤_ max is the energy budget. Since the optimum is attained when one of the constraints is active, the optimization program can be solved by a line search over μ (or τ) over the constraints curves, provided that r can be estimated at any given point. In our simulations we compute the function r for specific pairs of values (μ,τ) and take the minimum over these pairs. The two terms in the objective function of the static optimization problem show the tradeoff on the choice of the intermediate target isostable (which is to be reached after a time τ). For instance, choosing an isostable close to the equilibrium can lead to a large pulse duration τ (second term), but a small convergence time of the free motion (first term). Furthermore, the function r also allows to understand the tradeoff between the energy spent and the convergence time, which is not straightforward using standard optimal control theory. To summarize, we derived a static optimization problem, which has the same solution as the dynamic optimization problem (Problem 1). In order to compute the solution, one needs an efficient computational procedure for evaluating r at a given point. 


 §.§ Computation of the Pulse Control Function

The eigenfunction s_1 can be estimated through Laplace averages (<ref>) and the function r is subsequently obtained since it is the composition of the eigenfunction s_1 with the flow. In particular, we can derive the following formula:

    r(x,μ,τ)  = lim_t̅→∞1/t̅∫_0^t̅ g∘ϕ(t, ϕ(τ,x, μ),0) e^-λ_1 td t 
    
    = lim_t̅→∞1/t̅∫_τ^t̅ g∘ϕ(t, x, μ h(·,τ)) e^-λ_1 (t-τ)d t ,

where λ_1 is the dominant Koopman eigenvalue, g∈ C^1 satisfies g(x^∗)=0, v_1^T ∇ g(x^∗) ≠ 0, v_1 is the right eigenvector of J(x^∗) corresponding to λ_1 and h(t,τ) is the step function defined in (<ref>). In practice, we choose g(x)= w_1^T (x- x^∗), where w_1 is the dominant left eigenvector of J(x^∗). Since λ_1 is real according to Assumption A3 and Proposition <ref>), we have 

    r(x, μ, τ)      = lim_t̅→∞ w_1^T(ϕ(t̅, x, μ h(·,τ)) - x^∗) e^-λ_1 (t̅-τ)
       ≈ w_1^T(ϕ(t̅, x, μ h(·,τ)) - x^∗) e^-λ_1 (t̅-τ)

where the time t̅ should be chosen large enough.  In this case the tolerance of the differential equation solver should be set to O(e^λ_1 (t̅-τ)). When only observed data are available, the eigenfunction — and therefore the function r — can be computed through dynamic mode decomposition methods (cf. <cit.>). This idea is illustrated in Appendix <ref>.



 §.§ Is the Control Space Rich Enough?
 
Throughout the paper, we assume that the problem has a solution in the form of a temporal pulse u(t) = μ h(t, τ). We will argue that in the case of monotone systems, this is not a restrictive assumption. First of all, if the system is globally asymptotically stable, then clearly we can converge to x^∗ by using a temporal pulse.

Assume now that the system is monotone with two exponentially stable equilibria x^∗ and x^∙ and basins of attraction (x^∗) and (x^∙), respectively. Let the system be defined on a forward-invariant set = (x^∗)⋃(x^∙). Assume also that x^∗≫ x^∙, which is typically fulfilled in many bistable monotone systems. If x_0∈(x^∗), then we can choose u = 0, which is a temporal pulse with τ =0. Consider now the case x_0 = x^∙. If there exists a control signal u^1∈_∞ driving the system from x^∙ to x^∗, then we have ϕ(t, x^∙, u^1) ≼ϕ(t, x^∙, μ), where u^1(t)≤μ for (almost) all t. At a time τ, the flow ϕ(τ, x^∙, u^1) will be in the vicinity of x^∗ and in the basin of attraction of (x^∗). The flow ϕ(τ, x^∙, μ) will also be in the basin of attraction of x^∗. Indeed, if ϕ(τ, x^∙, μ) ∈(x^∙), then [x^∙, ϕ(τ, x^∙, μ)]∈(x^∙), which contradicts that ϕ(τ, x^∙, u^1) ∈(x^∗) and ϕ(τ, x^∙, u^1) ∈[x^∙, ϕ(τ, x^∙, μ)] (cf. <cit.>). Hence if we can switch from x^∙ to x^∗ with a control signal u(t), then we can switch with a temporal pulse. Finally the case x_0 ∈(x^∙) is treated in a similar manner by first allowing the trajectory to converge to a neighborhood x^∙ with u^2 =0 and then applying the argument above.

This discussion shows that using temporal pulses in the case of monotone systems does not restrict the space of feasible problems. However, we can strengthen the argument by showing that constant controls are optimal in the absence of energy constraints in Appendix <ref>. 




§ PROOF OF THE MAIN RESULT
 
In order to prove the main result, we would need to use the properties of the pulse control function r, which we present in the following lemma.


Let the system (<ref>) satisfy Assumptions A1–A5. Then r is a C^1 function on its effective domain (r). Furthermore, for all (x,μ,τ)∈(r)

(i)  ∂_x    r(x,μ,τ) ≫ 0,   ∂_μ r(x,μ,τ) > 0  and  ∂_τ r(x,μ,τ) > λ_1 r(x,μ,τ)

(ii) If r(x,μ,τ)≤ 0, then ∂_τ r(x,μ,τ)> 0;

(iii) If f(x, ν)≽ 0, then ∂_τ r(x,μ,τ)>0 for all finite τ> 0 and μ > ν.
 

(o) First, we show that under the assumptions above for all  t >0 and μ_1 > μ_2, we have

    ϕ(t, x, μ_1) ≻ϕ(t, x, μ_2), 
    	s_1(ϕ(t, x, μ_1)) > s_1(ϕ(t, x, μ_2)).

Due to monotonicity, we have that ϕ(t,x, μ_1) ≽ϕ(t,x, μ_2) for all  t >0  and  μ_1 > μ_2. All we need to show is that ϕ(t,x, μ_1) ϕ(t,x, μ_2) for all finite t>0. At t=0, the time derivatives of the flow are equal to f(x,μ_1) and f(x,μ_2). Since f(x,μ_1) ≻ f(x,μ_2) (Assumption A5), there exists a T>0 such that ϕ(t,x, μ_1) ≻ϕ(t,x, μ_2) for all t<T. If for some T we have that ϕ(T,x, μ_1) = ϕ(T,x, μ_2) and ϕ(t,x, μ_1) ≻ϕ(t, x, μ_2) for all t <T, then for some index i we have

    d ϕ_i(t,x, μ_1)/d t|_t =T <  d ϕ_i(t,x, μ_2)/d t|_t =T.

This implies that f_i(ϕ(T, x, μ_1),μ_1) < f_i(ϕ(T, x, μ_2), μ_2), which together with ϕ(T, x, μ_1) = ϕ(T,x, μ_2) contradicts Assumption A5. Therefore, ϕ(t,x, μ_1) ≻ϕ(t,x, μ_2) for all finite t> 0. Due to Assumption A4 we have that ∇ s_1(x )≫ 0, which in particular means that s_1(x) > s_1(y) for all x≻ y, and (<ref>) follows from (<ref>).

(i)  The flow is continuously-differentiable for constant control signals since f(x, u)∈ C^2 (Assumption A1), and hence r(x,μ,τ) = s_1(ϕ(τ,x,μ)) is a C^1 function. 

For x ≻ y, where  (x,μ,τ) and (y,μ,τ)∈(r), we have s_1(ϕ(τ, x,μ)) > s_1(ϕ(τ, y, μ)) due to monotonicity (Assumption A3) and Assumption A4. Hence ∇_x r(x, μ,τ)≫ 0.

For μ > ν, where  (x,μ,τ) and (x, ν,τ)∈(r), we have s_1(ϕ(τ, x, μ)) > s_1(ϕ(τ, x, ν)) due to monotonicity (Assumption A3) and point (o). Hence ∂_μ r(x, μ,τ)> 0. 

Finally,  ∂_τ r(x, μ,τ)> λ_1 r(x,μ, τ) follows from:

    ∂_τ r(x,μ, τ)=d s_1(ϕ(t,x,μ))dt|_t =τ = 
    ∇ s_1(ϕ( τ, x,μ))^T f(ϕ(τ, x,μ), μ) > 
    ∇ s_1(ϕ( τ, x,μ))^T f(ϕ( τ, x,μ), 0) = 
    λ_1 s_1(ϕ(τ, x, μ))=λ_1 r(x,μ, τ),

where the inequality is due to Assumption A4 and A5, and the following equality is due to (<ref>). 

(ii) This follows directly from point (i). 

(iii) This proof employs a fairly standard technique in monotone system theory. First note that f(x,ν)≽ 0 and Assumption A5 imply that f(x,μ)≻ 0. Consider a perturbed system ż = f(z,μ) +1/n, where  is a vector of ones and n is a positive integer. Let the flow of this system be ϕ_n(t,x,μ). We have that f(x,μ) +1/n ≫ f(x,μ)≻ 0. 	Since this is the derivative of the flow with respect to time around t=0, then ϕ_n(δ, x, μ)≫ x for a sufficiently small positive δ. Now for t > δ, we have

    ϕ_n(t, x, μ)  = ϕ_n(t-δ, ϕ_n(δ, x, μ), μ) ≫ϕ_n(t-δ, x, μ),
		
and hence ϕ_n(t, x, μ)  ≫ϕ_n(ξ, x, μ)	for any t>ξ. This conclusion holds for all n>0. With n→∞, we have that ϕ_n(t,x,μ)→ϕ(t,x,μ). Therefore, ϕ(t,x, μ) ≽ϕ(ξ, x, μ) for all finite t>ξ and ∂_τϕ(τ, x, μ) ≽ 0.  If the equality ∂_τϕ(τ, x, μ) = 0 is attained, then ϕ(τ, x, μ) is an equilibrium of the system ẋ = f(x, μ). This is impossible since an exponentially stable equilibrium cannot be reached in finite time τ due to uniqueness of solutions. Hence we have ϕ(t,x, μ) ≻ϕ(ξ, x, μ) for all finite t>ξ and, using Assumption A4, we obtain s_1(ϕ(t,x,μ)) > s_1(ϕ(ξ, x, μ)) for all finite t>ξ. The result follows.


If Assumptions A4 and A5 do not hold, then all the inequalities in Lemma <ref> are not strict. For instance, we have that  ∂_x r(x,μ,τ) ≻ 0, ∂_μ r(x,μ,τ) ≥ 0  and ∂_τ r(x,μ,τ) ≥λ_1 r(x,μ,τ) in point (i). We present additional properties of the function r in Appendix <ref>. Now we can present the proof of our main result.

Theorem <ref>
It is straightforward to verify that all the constraints and optimization objective are the same for Problem 1 and problem (<ref>). Hence, by construction the first part of the statement is fulfilled.

According to the constraint (<ref>), we have that r(x_0,μ,τ)<0, which implies the following chain of inequalities

    ∂_τ (ln(|r(x_0,μ,τ)| e^|λ_1| τ)) = 
    ∂_τ (|r(x_0,μ,τ)| e^|λ_1| τ)/|r(x_0,μ,τ)| e^|λ_1| τ = 
    -∂_τ (r(x_0,μ,τ)) ·  e^|λ_1| τ +  |λ_1|  |r(x_0,μ,τ)| e^|λ_1| τ/|r(x_0,μ,τ)| e^|λ_1| τ < 
    λ_1  |r(x_0,μ,τ)| e^|λ_1| τ +  |λ_1|  |r(x_0,μ,τ)| e^|λ_1| τ/|r(x_0,μ,τ)| e^|λ_1| τ = 0

where the inequality follows from Lemma <ref>. Hence, the derivative of the objective function in (<ref>) with respect to τ is negative. Finally, ∂_μln(|r(x_0,μ,τ)| e^|λ_1| τ) is also negative according to Lemma <ref>. Hence, if there is a feasible point, the constraints (<ref>), (<ref>) are reached in order to minimize the objective. 






§ APPLICATIONS



 §.§ Implementation Details

We implemented our computational procedures both in python (using lsoda ordinary differential equation (ODE) solver) and Matlab (using ode15s ODE solver). We run our computational algorithm on a laptop equipped with a 4 core Intel i7 processor running at 2.4 GHz and 8 GB of RAM, however, we did not explicitly parallelize the computations. Computing one value of the function r is equivalent to computing one trajectory of the system, albeit with high precision (we set relative tolerance of the solvers to e^-10 – e^-14). 



 §.§ Closed-Loop Switching in Generalized Repressilator

The eight species generalized repressilator is an academic example, where each of the species represses another species in a ring topology (cf. <cit.>). The corresponding dynamic equations for a symmetric generalized repressilator are as follows: 

    ẋ_1    = p_1 1^0/1 + (x_8/p_1 2^0)^p_1 3^0 + p_1 4^0 - p_1 5^0 x_1 + u, 
    ẋ_i    = p_i 1^0/1 + (x_i-1/p_i 2^0)^p_i 3^0 + p_i 4^0 - p_i 5^0 x_i, ∀ i = 2,… 8,

where p_i 1^0 = 40, p_i 2^0 = 1, p_i 3^0 = 2, p_i 4^0 = 1, and p_i 5^0 = 1. This system has two exponentially stable equilibria x^∗ and x^∙ and is monotone with respect to the cones _x = P_x ^8 and _u =, where P_x= ([1, -1, 1, -1, 1, -1, 1, -1]). We have also x^∙≼__x x^∗. It can be shown that the unforced system is strongly monotone in the interior of ^8 for all positive parameter values. We consider here the problem of switching the system from one equilibrium x^∙ to the other equilibrium x^∗ and we can verify that there exist pulse control signals that induce such a switch.

In the left panel of Figure <ref>, we plot the level set of the function _ conv(x^∙, μ,τ, 10^-2) = 1/|λ_1|ln(r(x^∙, μ, τ)/10^-2) + τ, the level set _ max = 100 of the function μτ = _ max, and the level set r(x^∙,μ,τ) = 0. The last two are related to the constraints of the static optimization program (<ref>). Note that r is computed with the dominant eigenfunction associated with the target equilibrium x^∗.  We also note that the function _ conv can escape to -∞ around the level set r(x^∙, μ,τ) ≈ 0. This is not a conflict with the interpretation of  the function _ conv, since it represents the convergence time only if the value of |r(x^∙,μ,τ)| is larger than 10^-2. Otherwise the term 1/|λ_1|ln(|r(x^∙, μ, τ)|/10^-2) is negative, and the computational results are meaningless. This also explains why the level sets of _ conv appear to have the same asymptotics as the level set r(x^∙, μ,τ) ≈ 0 in Figure <ref>. 

Our goal is to compare the open-loop (proposed by <cit.>) and closed-loop solutions to the switching problem subject to perturbations of parameters p_i j^0. We consider two settings for the simulation. In both settings, we compute the control signals based on the nominal model (<ref>) with the parameter values p_i j^0, but the simulations are obtained with two sets of (exact) parameter values:

Setting A.  We set p_i 1^A = 50 for odd i.


Setting B. We set p_i 1^B = 30 for odd i. 

The Euclidean distance between the nominal initial point and the actual initial point in Setting A and B is equal to 0.025 and 0.031, respectively. In order to compute an open-loop optimal control policy based on the nominal model (i.e. with parameter values p_i j^0), one can solve the static optimization program (<ref>). The plots in Figure <ref> also offer a graphical solution to the problem and a depiction of possible tradeoffs in the problem. In our case, the optimal solution lies at the intersection of the constraint curves (i.e. energy budget curve and level set r(x^∙,μ,τ) = 0). 

In our simulations, we pick a pair (μ^0, τ^0) lying near the zero level set of r below the level set _ max = 100. This is not an optimal solution for the energy budget _ max = 100, however, we pick a solution with a lower energy expenditure and a larger time τ (and hence larger convergence time) in order to have a possibility to react to the obtained measurements in the case of the closed-loop setting. We take τ^0 = 20 and compute μ^0 minimizing the time _ conv on a uniform grid of 100 points in [2, 10], which gives the value μ^0 = 3.53.

For the closed-loop control, we take the same initial pair (τ^0, μ^0).  In both setting A and B, we update the control signals every t_ samp = 2. For each update, we decrease the time of the pulse by t_ samp and we decrease the available energy budget by subtracting the energy already consumed. We then compute the values of the function r with a fixed τ = 20 - N t_ samp, where N is the number of previous updates and we choose the value μ on a uniform grid of 100 points in [2, 10], which minimizes _ conv with r(x^∙, μ, τ) < 0. 

The simulation results are depicted in the center and right panels of Figure (<ref>). In Simulation A, the system converges to the target equilibrium faster than the nominal one (i.e, with parameters p_i j^0) and the closed-loop solution saves energy and limits the overshoot in comparison with the open-loop solution. In Simulation B, the opposite occurs and all the energy budget is spent. In this case, the closed-loop solution allows the switch, while the open-loop (i.e. <cit.>) does not. 



 §.§ Synchronization of Cardiac Cells


Besides the switching problem, our approach can be used for more general problems. For example, a problem of importance in biology is to synchronize an ensemble of systems (e.g. cells) with a common control input. In this subsection, we propose to use the function r in order to compute a train of temporal pulses and solve the synchronization problem.  Provided that the function r and the state x^j of each system are known, we design a pulse that drives every system toward isostables. The pair (μ^∗,τ^∗) is chosen in such a way that the corresponding pulse minimizes the maximum time delay between the different systems, and it follows from (<ref>) that this pair can be computed as

    (μ^∗,τ^∗) = _(μ,τ) ∈^2ln(max_j r(x_j,μ,τ)/min_j r(x_j,μ,τ)).

From a practical point of view, we note that the values of the function r might be known only for some (x,μ,τ) (especially if it is computed from data). In this case, the value r at the different states x_j of the systems can be computed by interpolation, and the optimal pair (μ^∗,τ^∗) will be picked among the pairs (μ,τ) for which the value of r is known. 



We apply this closed-loop control (see Algorithm <ref>) to synchronize FitzHugh-Nagumo systems (cf. <cit.>), which have been proposed as simple models of excitable cardiac cells. This example is motivated by the synchronization of cardiac cells (i.e. defibrillation) and is directly inspired by the study by <cit.> proposing optimal defibrillation strategies. We consider here 100 FitzHugh-Nagumo cells described by the dynamics (see also <cit.>)

    V̇    = 0.26 V(V-0.13)(1-V)-0.1 V w 
    ẇ    = 0.013 (V-w)

where V is the membrane potential and w is a recovery (gating) variable. The values of r were computed a priori on a 20 × 20 grid for (V,w) ∈[0,2]× [0,2] and on a 51 × 41 grid for (μ,τ) ∈ [0,0.5] × [10,50]. The time between two successive pulses is T_p=70 and the initial conditions are randomly distributed on [0,2]×[0,2]. The maximum time delay between the cells (computed with (<ref>)) after each pulse is shown in Figure <ref>(a) for the input obtained with the closed-loop control (optimal pairs (μ^∗,τ^∗)) and for periodic pulse trains. The best performance is obtained with the closed-loop control. As shown in Figure <ref>(b), the optimal pairs (μ^∗,τ^∗) are not identical at each iteration, since they depend on the states of the cells, which motivates the use of closed-loop control. The first pulses correspond to a maximum value τ^*=50, but μ^* takes intermediate values in the interval [0,0.5]. In particular, small values of μ are needed to obtain a fast convergence rate. We observe in Figure <ref>(a) that periodic pulse trains also synchronize the cells, but slower than our closed-loop approach. A periodic pulse train with maximum values (μ,τ)=(0.5,50) yields a slow rate of convergence (red curve), while smaller values (μ,τ) yield very large delays for the first iterations (green curve). Clearly, the optimal approach by <cit.> outperforms our method in terms of convergence and time delays, however, our optimal control policy is easier to implement. Furthermore, we can parametrize our control signal with different (non-constant) basis functions.



The proposed closed-loop control could be easily adapted to incorporate additional constraints (e.g. maximum energy) and provides a solution to the synchronization problem. This solution is a compromise between a simple but not optimal periodic pulse train and the complex exact solution of the optimal control problem established by <cit.>. Future work could extend these preliminary results to more realistic cases, for instance where not all the states are observable.



§ CONCLUSION


In this paper, we studied a switching/convergence problem for monotone systems. Our solution reduces a dynamic optimization problem to the computation of the time-independent function r, which is defined using the Koopman operator. The properties of the function r lead to straightforward solutions to a tradeoff between the convergence time and the energy budget. This approach can potentially be extended beyond monotone systems and switching/convergence problems. In this paper, we illustrate the possible benefits of a closed-loop solution for the switching problem. We also apply our framework to the synchronization of cardiac cells represented by non-monotone FitzHugh-Nagumo models. In this paper, we have not addressed partial state observability and/or partial controllability issues. This constitutes one of the future work directions.






§ CONSTANT CONTROL SIGNALS ARE OPTIMAL FOR A MINIMUM-TIME PROBLEM
 
Consider the following optimal control problem over bounded measurable control signals: 

    V(z, μ, β)    = inf_[ τ, u ∈_∞([0,μ]) ]τ, 
     subject to    (<ref>),  x(0) = z, 
       x(τ) ∈_β = {y∈^n | s_1(y) = β},

where  s_1(x) is a C^1 increasing dominant eigenfunction defined on the basin of attraction of x^∗. Under our assumptions, the solution to this problem is surprisingly straightforward.

	Let the system (<ref>) satisfy Assumptions A1 – A3. Then 
	
	(i) If s_1(z) < β, then the optimal solution to (<ref>), if it exists, is u(t) = μ for all t∈[0, τ];
	
	(ii) If s_1(z) ≥β, then the optimal solution to (<ref>) is u(t) = 0 for all t≥0. 	
	 
	
		(i) Let u^0(t) = μ for all t>0, and u^δ(t) be any admissible control signal, then u^0(t) ≽ u^δ(t) for all t∈[0, τ]. Then by monotonicity we have ϕ(t, z, u^0(·))≽ϕ(t, z, u^δ(·)), which leads to s_1(ϕ(t, z, u^0(·)))≥ s_1(ϕ(t, z, u^δ(·))) for all t≥ 0 due to Proposition <ref>.  Hence
		
    β   = s_1(ϕ(τ, z, u^0(·)))≥  s_1(ϕ(τ, z, u^δ(·))),
    β   > s_1(ϕ(t, z, u^0(·))) ≥  s_1(ϕ(t, z, u^δ(·)))  for  t<τ,

		which implies that the target set _β^- is reached with u^0(·) at least as fast as with any other admissible control signal u^δ(·). Therefore, the control signal u^0(t) = μ is an optimal solution of the problem.
		
		(ii) The proof is similar to the point (i). 
		
This result justifies our use of temporal pulses to solve convergence problems for monotone systems. The problem (<ref>) has a direct relation to the function r. In particular, s_1(x)<β implies that r(x, μ, V(x, μ, β)) = β provided that the problem (<ref>) has a solution. On the other hand, there might be some values τ V(x, μ, β) such that  r(x, μ, τ) = β. In general, if s_1(x)<β, then 

    V(x, μ, β) = min{τ∈ | r(x,μ,τ) = β},

provided that the solution exists. 
In particular, if the premise of the points (ii) and (iii) in Lemma <ref> holds then V(x, μ, β) = τ if and only if r(x,μ,τ) = β, which again justifies our use of the pulse control function in these cases.


§ LEVEL SETS OF THE FUNCTION R
 
If the premise of the points (ii) and (iii) in Lemma <ref> hold, then the level sets of r are the graphs of strictly decreasing functions, a result which simplifies their computations as discussed in <cit.>. The algorithms describing this procedure were developed by <cit.> and <cit.>.

 
	Let the system (<ref>) satisfy Assumptions A1, A4, A6, A7, and α > 0. 

	(i) If s_1(x) < -α, then {(μ,τ) ∈^2| r(x, μ,τ) = -α} is a graph of a strictly decreasing function, i.e., this set does not contain pairs (μ_1, τ_1) (μ_2, τ_2) such that μ_1≤μ_2 and τ_1 ≤τ_2;

	(ii) If s_1(x) < α and f(x,η)≽ 0, then {(μ,τ) ∈^2 | r(x, μ,τ) = α} is a graph of a strictly decreasing function for μ>η.
	
	(i) First, we make sure that the function r(x,·,·) is increasing.  This is guaranteed if ∂_μ r(x,μ,τ), ∂_τ r(x,μ,τ) are positive. 
		
		Due to Lemma <ref> we have that ∂_τ r(x, μ, τ) >0 as long as r(x, μ, τ)≤ 0, and that ∂_μ r(x,μ,τ)>0.
		Now since ∂_τ r(x, μ, τ) >0 (for all (μ,τ) such that r(x,μ,τ)=-α<0), the implicit function theorem implies that there exists a function τ = g(μ, α) such that r(x, μ, g(μ, α)) = -α. Furthermore, the function g is C^1 in μ and ∂_μ g = - ∂_μ r(x, μ, τ) / ∂_τ r(x, μ, τ) < 0 in the neighborhood of the level set { (μ, τ)  | r(x, μ, τ) =-α}. Therefore the level set {τ, μ | r(x, μ, τ) =-α} is a graph of a strictly decreasing function in μ. It also directly follows that  the level set { (μ , τ) | r(x, μ, τ) =-α} is a graph of a strictly decreasing function in τ.
		
		(ii) We have that ∂_μ r(x,μ,τ) > 0. Now positivity of ∂_τ r(x,μ,τ) for positive r(x,μ,τ) follows from the point (iii) in Lemma <ref>. The rest of the proof follows on similar lines as in (i).
		


§ PURELY DATA-BASED CONTROL SETTING

One can compute r from sampled data points (e.g. obtained from experiments). If data points
 
    z[k]=ϕ(τ+k T_s, x, μ h(·, τ)) 
    = ϕ(k T_s,ϕ(τ, x, μ),0)

 are given (with the sampling time T_s), then the DMD algorithm can be used to compute s_1(ϕ(τ, x, μ ))=r(x,μ,τ). We illustrate this in the following subsections.


 §.§ Dynamic mode decomposition algorithm
 
The dynamic mode decomposition algorithm can be used to estimate the eigenfunctions of the Koopman operator from data. Assume that N snapshots of m trajectories of the (unforced) system are given, i.e.

    z^j[k] = g(ϕ(k T_s,x_0^j,0)) ,   j=1,…,m , k=1,…,N

where T_S is the sampling time, x_0^j ∈ℝ^n is the initial condition, and g is an observable (typically measuring one of the states). The algorithm yields the so-called DMD modes and eigenvalues, which are related to the eigenfunctions and eigenvalues of the Koopman operator. It is described in Algorithm <ref>. We will not elaborate further on the details of the method and refer the reader to <cit.>. We only note that the matrix Y X^†, where X^† is the pseudo-inverse of X, is a finite-dimensional approximation of the Koopman semigroup. Therefore the eigenvectors of Y X^† can be used to estimate the eigenfunctions of the Koopman operator.




 §.§ Application of the DMD Algorithm to Toxin-antitoxin System
 

Consider the toxin-antitoxin system studied by <cit.>. We will use the model and parameter values as described in <cit.>. 

The system is bistable with two exponentially stable equilibria x^∙ and x^∗, but not monotone with respect to any orthant. However, it was established by <cit.> that it is eventually monotone, i.e. the flow satisfies the monotonicity property after some initial transient. The level sets of the function r(x^∙, μ, τ) (in the space of the parameters (μ, τ)) were computed by <cit.> and it appears that these sets are monotone curves although the system is not monotone. 

Since we cannot guarantee that the function r is increasing, we need to develop a different computational approach applicable to a broader class of systems. One possibility is to compute the function r from data. To do so, we used the DMD algorithm with 25 time series of the form (<ref>) (each of which corresponds to a different pair (μ,τ) ∈ [24,50] × [3,25]). For each time series, we used only 4 snapshots over the time interval [τ+12.5, τ+50] (i.e. k ∈{1,2,3,4} with T_s=12.5 in (<ref>)). The level sets of r computed with the DMD algorithm are accurate and similar to the ones obtained with Laplace averages (Figure <ref>). We stress that this example is only a proof-of-concept, and further research is required to automate the application of DMD to this problem.

