Peak Estimation of Rational Systems using Convex Optimization

This paper presents algorithms that upper-bound the peak value of a state function along trajectories of a continuous-time system with rational dynamics. The finite-dimensional but nonconvex peak estimation problem is cast as a convex infinite-dimensional linear program in occupation measures. This infinite-dimensional program is then truncated into finite-dimensions using the moment-Sum-of-Squares (SOS) hierarchy of semidefinite programs. Prior work on treating rational dynamics using the moment-SOS approach involves clearing dynamics to common denominators or adding lifting variables to handle reciprocal terms under new equality constraints. Our solution method uses a sum-of-rational method based on absolute continuity of measures. The Moment-SOS truncations of our program possess lower computational complexity and (empirically demonstrated) higher accuracy of upper bounds on example systems as compared to prior approaches.


I. INTRODUCTION
Peak estimation is the practice of finding extreme values of a state function p along trajectories x(t) of a dynamical system that evolve starting from an initial set X 0 .Instances of peak estimation (extremizing p(x(t))) include finding the maximum speed of an aircraft, height of a rocket, concentration of a chemical, and current along a transmission line.This work focuses on peak estimation in the case of rational continuous-time dynamics for a state x ∈ R n where: where the rational dynamics f can be represented as The expression in (1) is a sum-of-rational dynamical system in terms of polynomials f 0 , N ℓ , and D ℓ (with L finite).Applications of peak estimation for rational systems include systems include finding maximal concentrations in chemical reaction networks with Michaelis-Menten kinetics or yeast glycolysis, velocities in rigid body kinematics (manipulator equation with rational friction models), and occupancies in network queuing models [1], [2].Refer to [1] for a detailed survey of applications of rational systems, as well as a formulation of algebraic analysis techniques to establish system properties such as parameter identifiability and controllability.
The rational-dynamics peak estimation task considered in this work (maximizing a state function p along system trajectories x(t | x 0 ) evolving in a state space X ∈ R n starting from X 0 ⊆ X with a time horizon of [0, T ]) is described in Problem 1.
All of the previously mentioned applications of LP in dynamical systems analysis and control (in the context of continuous state spaces) use a Lipschitz assumption in dynamics in order to prove that there is no relaxation gap between the infinite-dimensional LP and the original finitedimensional nonconvex program.The rational dynamics in (1) may fail to be globally Lipschitz (within the domain of [0, T ] × X), and therefore falls into the theory of nonsmooth dynamical systems [23], [24].This work utilizes a sumof-rational representation from [25] in order to cast (3) as an LP in measures, and uses the theory of nonsmooth Liouville equations from [23] to prove equivalence of optima under compactness, trajectory-uniqueness, and positivity assumptions.Prior work in SOS-based analysis of rational functions includes clearing to common denominators [26] (under positivity), and adding new variables to represent the graph of rational functions as equality constraints [27], [28].
The main contributions of this work are: • A measure LP formulation for peak estimation of rational dynamical systems based on the theory of [25] (for sum-of-rational static optimization).• A proof that there is no relaxation gap between the finite-dimensional nonconvex problem (3) and the ob-jectives of infinite-dimensional LPs • Quantification of the computational complexity in terms of sizes of the Positive Semidefinite (PSD) matrices in SDPs.
• Experiments demonstrating the upper-bounding of the true peak on rational dynamical systems.This paper has the following structure: Section II reviews preliminaries such as notation and occupation measures.Section III formulates a sum-of-rational-based measure LP for the peak estimation problem in (3).Section IV introduces and applies the moment-SOS hierarchy to obtain a nonincreasing sequence of upper-bounds to the true peak value P * .Section V performs peak estimation on example rational dynamical systems.Section VI concludes the paper.
An extended version of this conference paper is available at [29].The extended version includes a proof of strong duality, as well as details about SOS programs for the preexisting rational-peak-estimation methods [26], [27].

A. Notation
The set of n-dimensional indices with sum less than or equal to a value d is The set of polynomials with indeterminate x is R[x], and the subset of polynomials with degree at most The set of continuous (continuous and nonnegative) functions over S is C(S) (C + (S)).The set of signed (nonnegative) Borel measures over a set S is M(S) (M + (S)).The measure of a set A ⊆ S w.r.t.µ ∈ M + (S) is µ(A).The sets C + (S) and M + (S) possess a bilinear pairing ⟨•, •⟩ that acts by Lebesgue integration: g ∈ C + (S), µ ∈ M + (S) : ⟨g, µ⟩ = S g(s)dµ(s).This bilinear pairing is an inner product between C + (S) and M + (S) when S is compact (in which M + (S) can be canonically identified as the dual of C + (S)), and the pairing can be extended to integration between C(S) (and sets of more general measurable functions) and M(S).Given two measures µ 1 ∈ M + (S 1 ), µ 2 ∈ M + (S 2 ), the product measure µ 1 ⊗ µ 2 is the unique measure satisfying and µ is a probability measure if this mass is 1.The Dirac delta supported at a point s ′ (δ s=s ′ ) is the unique probability measure such that ∀g ∈ C(S) : ⟨g, δ s=s ′ ⟩ = g(s ′ ).The adjoint of a linear operator L is L † .

B. Occupation Measures
An occupation measure is a nonnegative Borel measure that contains all possible information about the behavior of (a set of) trajectories of a given dynamical system.
For a given initial condition x 0 ∈ X 0 , the occupation measure µ The (t * , x * 0 )-occupation measure µ x(•) in ( 4) can also be understood in terms of its pairing with arbitrary continuous (measurable) functions: (5) Occupation measures µ x(•) in ( 4) may be defined over a distribution of initial conditions µ 0 ∈ M + (X 0 ) (with x 0 ∼ µ 0 ): Any trajectory of (3b) satisfies the conservation law, The conservation law in ( 7) is a Liouville equation, and can be written in terms of the initial measure The ∀v imposition in equation ( 8) can be written equivalently in shorthand form (with L † f as the adjoint linear operator of L) as Any triple (µ 0 , µ p , µ) that satisfies ( 9) is a relaxed occupation measure; the class of relaxed occupation measures may be larger than the set of superpositions (distributions) of occupation measures arising from trajectories.

III. RATIONAL LINEAR PROGRAM
This section will present convex infinite-dimensional LP to perform peak estimation of rational systems.

A. Assumptions
We will begin with the following assumption: A1: If a trajectory satisfies x(t | x 0 ) ̸ ∈ X for some x 0 ∈ X 0 , then x(t ′ | x 0 ) ̸ ∈ X for all t ′ ≥ t.Further assumptions will be added as needed.
Remark 2: Assumption A1 is a non-return assumption in the style of [30].
Proof: By Theorem 3.1 of [24], imposition of assumption A4 ensures that every relaxed occupation measure (µ 0 , µ p , µ) is supported on the graph of (a superposition of) trajectories of (3b).Compactness (A2) and (lower semi-) continuity (A3) are necessary to invoke arguments used by Theorem 2.1 of [3], using the non-smooth Theorem 3.1 of [24] rather than a Lipschitz assumption on dynamics.

C. Absolute Continuity Formulation
We will use the sum-of-rationals framework of [25] in order to express (10) in a form more amenable to numerical computation, using the moment-SOS hierarchy of SDPs.This sum-of-rationals framework uses the notion of absolute continuity of measures.
We now impose the following assumption: A5: Each function D ℓ is strictly positive over [0, T ] × X.
Proposition 9: The measures ν ℓ have finite densities The Lie derivative in (10b) can be expanded using ( 11) as (∀v ∈ C 1 ([0, T ] × X) with assumption A5 in place) Problem 10 uses the sum-of-rational framework to pose a measure peak estimation LP: Problem 10: Find an initial measure µ 0 , a relaxed occupation measure µ, a peak measure µ p , and a set of per-rational measures {ν ℓ } L ℓ=1 to supremize (14f) Corollary 11: Under A1-A5, the objective values of Problem 1 and 10 will be equal with P * = p * r .Proof: Theorem 5 ensures that P * = p * under A1-A4.The finite nature of the densities 1/D ℓ from Proposition 9 (proved by Theorem 2.1 of [25]) ensures that the absolutecontinuity-based construction process in (11) will result in each ν ℓ having identical support to µ.As a result, (µ 0 , µ p , µ) from (14b) will form a relaxed occupation measure, thus ensuring no relaxation gap by Theorem 3.1 of [24] (used in Theorem 5).
Remark 12: If A5 is not imposed, then it is possible for the measure ν ℓ from (11) to be unconstrained (with D ℓ = 0 at some (t, x)).It therefore cannot be presumed that the supports of ν ℓ and µ are identical.These degrees of freedom in ν ℓ could allow for the strict (and possibly unbounded) upper-bound p * r > p * .

IV. FINITE-DIMENSIONAL TRUNCATION
Program (15) must be discretized into a finite-dimensional program in order to admit tractable numerical solutions by computational means.We will first introduce the moment-SOS hierarchy of SDPs, and then use this hierarchy in order to perform finite-dimensional truncations of (15).

A. Sum of Squares Background
The set of all SOS polynomials in indeterminates x is marked as Σ[x] ⊂ R[x], and the boundeddegree subset of SOS polynomials with degree less than or equal to 2k is Σ[x] ≤2k .The set of SOS polynomials is a strict subset of the cone of nonnegative polynomials, with equality holding only in the cases of univariate polynomials, general quadratics, or bivariate quartics [31], [32].To each SOS polynomial θ, there exists a (nonunique) tuple of a size s ∈ N, a polynomial vector m(x) ∈ N s , and a PSD Gram matrix Q ∈ S s + such that θ(x) = m(x) T Qm(x).When x has dimension n and m(x) is chosen to be the vector of monomials of degrees 0 to k, the size s is n+d d .Testing membership of a polynomial in the SOS cone can therefore be done using SDPs [33].
A Basic Semialgebraic (BSA) set is a set defined by a finite number of bounded-degree polynomial inequality constraints.For any BSA set K = {x | g j (x) ≥ 0, j ∈ 1 . . .N g }, the Weighted Sum of Squares (WSOS) cone Σ[K] is the class of polynomials θ ∈ R[x] that admit the following representation in terms of SOS polynomials (σ 0 , {σ j }): Every compact set whose bounding radius R is known may be rendered ball-constrained by appending the redundant constraint R −∥x∥ 2  2 to the description of K. Every positive polynomial over a ball-constrained set K is also a member of Σ[K] (Putinar Positivestellensatz [34]).The multipliers σ j from (16) certifying this positivity may generically have degrees that are exponential in n and degree of θ [35].
The truncated WSOS cone Σ[K] ≤2k is the class of polynomials such that deg(σ 0 ) ≤ 2k and ∀j : deg(σ j ) ≤ 2k.The process of replacing a nonconvex polynomial inequality constraint with a WSOS constraint and increasing the degree until convergence is called the moment-SOS hierarchy.

B. SOS Program
We will impose the following constraints in order to utilize the moment-SOS hierarchy: A6: X 0 and X are ball-constrained BSA sets.
For a fixed degree k ∈ N and index ℓ ∈ 1 . . .L, let us define the following degree of: The degree-k SOS truncation of ( 15) is: Problem 16: Find a scalar γ and polynomials v, (18g) The following lemma ensuring finiteness of measure masses is required to prove convergence of Problem 16 to the optimal value of Problem 10 as k → ∞: Theorem 17: Under assumptions A1-A6, the finite truncations will converge from above as lim k→∞ d * r,k = P * Proof: We first note that P * = d * r under A1-A5 by Corollary 11 using Theorem 15 (strong duality) and Lemma 14 (finite mass).
After noting that the masses of all feasible measures are bounded by Lemma 14, convergence in objective to d * r = P * is proven by Corollary 8 of [36].

C. Computational Complexity
The dominant-size Gram PSD constraint of program (18) occurs at (18d), and has size n+1+2k+2⌊deg(f0)/2⌋ n+1 .When using an interior point method, the scaling of (18) therefore grows as O(n 6k ) (nominally) or exponentially (in a degenerate case from Proposition 6 of [37]).Further complexity reductions such as symmetry and term sparsity may be employed if present to reduce computation time, but exploiting sparsity may lead to different finite-degree SOS optimal values.
V. NUMERICAL EXAMPLES Julia code to generate all examples in this work is publicly available online 1 .SOS programs were posed using a Correlative-Term-Sparsity interface (CS-TSSOS) [38], [39].The SOS programs were converted to SDPs using JuMP [40].All SDPs were then solved by Mosek 10.1 [41].

A. Two-Species Chemical Reaction Network
The first example involves analysis of a chemical reaction network.The states (x 1 , x 2 ) represent the nonnegative concentrations of the two species.The species undergo degradation at a linear rate, and promotion according to Michaelis-Menten kinetics (therefore possessing a globally asymptotically stable equilibrium point in the nonnegative orthant [42]).The relevant dynamics evolving in X = [0, 1] 2 over a time horizon of T = 6 are: We note that the unique equilibrium point of ( 19) occurs at x eq = [0.3203,0.7027].Both of the denominators in (19) are positive over X, thus satisfying assumption A5.Peak estimation is performed to bound the maximal value of p(x) = x 2 for trajectories beginning in the disc initial set of A lower-bound on p(x) = x 2 acquired from gridded numerical ODE-sampling is 0.8157.Table I reports upperbounds acquired by finite-degree SOS truncations in (18), as well as bounds discovered by comparison methods [26], [27] detailed in Appendix II of [29] at the same degrees for v(t, x).Table II reports the time taken for Mosek to return solutions in Table I.All methods return an upper bound of P * ≤ 1 at degree k = 1.Our method (Problem 16) returns the lowest peak estimate at each degree k for this experiment.Degree k = 6 bounds (black dotted lines) and sample trajectories (colored curves) starting from X 0 are plotted in Figure 1.The k = 6 bounds from Table I are visually indistinguishable on Figure 1.The red dot marks the location of the stable equilibruim point x eq .

B. Three-state Rational Twist System
The second example involves a rational modification of the Twist system from Equation (37) of [30].The threestate rational twist system is expressed in terms of matrix parameters (A, B) to form the expression of (∀i ∈ 1..3): The polynomial Twist dynamics in [30] lacks the (0.5 + x 2 i ) denominators.Peak estimation of ( 21) occurs over the state space of X = [−1, 1] 3 in a time horizon of T = 6.The initial set X 0 is the two-dimensional box X | x3=0 .It is desired to   A lower bound on P * acquired through sampling (gridding the plane X 0 ) is P * > 0.3489.Sampled trajectories and the p(x) = x 2 3 level set at the k = 6 bound (Problem 16) for the rational Twist system are plotted in Figure 2. VI.CONCLUSION This paper presents a scheme to perform peak estimation of rational systems.The sum-of-rational optimization tech-nique from [25] is used to reduce the complexity of resulting moment-SOS-derived SDPs.This decomposition scheme is nonconservative if all denominator polynomials D ℓ are positive and a compact set is considered (generating assumptions A1-A6 in the dynamical systems setting).Effectiveness of our technique was demonstrated on example systems.
Future work includes investigating methods to reduce conservatism of assumptions A1-A6 (such as if D ℓ > 0 only along the graph of trajectories starting at X 0 ).Other options include applying our method towards rigid body kinematics, decomposing network structure in dynamics [43] through sparse sum-of-rational optimization (Theorem 3.2 of [25]), analyzing non-ODE rational system models, and applying the methods from [25] towards other dynamical systems problems [14]- [22].

Fig. 1 :
Fig. 1: Trajectories and k = 6 bounds for (19), along with position of the unique equilibrium point x eq upper-bound the peak value of p(x) = x 2 3 along trajectories of (21).Tables III and IV compile computed upper-bounds of p(x) along rational Twist trajectories and solution timing, in the same style as in Tables I and II.The Sum-of-Ratios program (18) returns a bound of P * ≤ 1 at degree k = 1.

TABLE II :
Timing (seconds) for TableI