Convex and Concave Relaxations for the Parametric Solutions of Semi-explicit Index-One Differential-Algebraic Equations

Abstract

A method is presented for computing convex and concave relaxations of the parametric solutions of nonlinear, semi-explicit, index-one differential-algebraic equations (DAEs). These relaxations are central to the development of a deterministic global optimization algorithm for problems with DAEs embedded. The proposed method uses relaxations of the DAE equations to derive an auxiliary system of DAEs, the solutions of which are proven to provide the desired relaxations. The entire procedure is fully automatable.

Appendices

Appendix A: Construction of u _ψ , o _ψ , u _f and o _f Relaxations

The state relaxations described in this article rely on the ability to compute functions u _f , o _f , u _ψ and o _ψ satisfying Definitions 4.8 and 4.6, as well as Assumption 4.1. A method for automatically generating and evaluating such functions based on McCormick’s relaxation technique [26] was described in [24], and was applied to a very similar application, namely relaxing the right-hand side functions of a system of ODEs, in Sect. 7.2 of that article. This section bridges the gap between [24] and the present article by providing precise definitions of the functions u _f , o _f , u _ψ and o _ψ through the algorithms defined in [24]. Accordingly, this section is not intended to be independent and should be read alongside [24]. Specifically, the notations $\mathcal{V}$, $\tilde{\mathcal{U}}$, $\tilde{\mathcal{O}}$ and Φ from [24] are used below.

Choose any i∈{1,…,n _x }. Then, letting $\mathcal{V}$ denote the collection of factors v ₁,…,v _m defining f _i (Definition 2.3) and letting the set Φ of Definition 15 in [24] be defined as Φ in Definition 4.3 of the present article, the complete algorithmic construction of the function u _f,i is specified by Definition 15 in [24] and the following definitions:

From Definition 15 in [24], the analogous definition of o _f,i is apparent. Given these definitions, the composite relaxation properties Definition 4.8 follow directly by construction and the rules of McCormick’s relaxations. The reader is referred to Theorems 8 and 14 in [24] for formal arguments.

The article [24] provides modest conditions (see Theorem 6) on the factorable representation of f (Definition 2.3), which guarantee continuity of u _f and o _f on $I\times P\times\mathbb{R}^{n_{x}}\times\mathbb {R}^{n_{x}}\times\mathbb{R}^{n_{y}}\times\mathbb{R}^{n_{y}}$, as well as the global Lipschitz condition of Assumption 4.1 (Condition 3). It is worth noting that these conditions do not imply a global Lipschitz condition on f, but they do imply a much less restrictive local one. Essentially, for fixed (t,p)∈I×P, the global Lipschitz condition on u _f and o _f is made possible by the intervals X(t) and Y(t). The construction of these functions involves mapping any arguments $(\mathbf{z}_{x}^{c},\mathbf{z}_{x}^{C},\mathbf{z}_{y}^{c},\mathbf {z}_{y}^{C})\in\mathbb {R}^{n_{x}}\times\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{y}}\times\mathbb {R}^{n_{y}}$ into X(t)×X(t)×Y(t)×Y(t) in a Lipschitz manner (using the $\operatorname{mid}$ function), so that Lipschitz continuity of u _f (t,p,⋅,⋅,⋅,⋅) and o _f (t,p,⋅,⋅,⋅,⋅) need only hold on this (compact) interval [24].

Now choose any i∈{1,…,n _y } and consider the functions u _ψ,i and o _ψ,i of Definition 4.6. The definition of these functions is slightly more complicated than that of u _f,i and o _f,i because ψ _i is only factorable with t, and hence C(t), constant. However, if the definition of ψ is modified, call it ψ′, so that $\hat{\mathbf{C}}\in\mathbb{R}^{n_{y}\times n_{y}}$ is taken as an independent augment and used in place of C(t), then $\psi_{i}'$ is clearly factorable and defined on

Let $\mathcal{V}$ denote the corresponding collection of factors v ₁,…,v _m and denote the elements of C by c _i , where $i\in\{1,\ldots,n_{y}^{2}\}$. Letting the set Φ of Definition 15 in [24] be defined as Φ′ above, the complete algorithmic construction of the function u _ψ,i is specified by Definition 15 in [24] and the following definitions:

Again, the analogous definition of o _f,i is apparent from Definition 15 in [24]. Given these definitions, the composite relaxation properties of Definition 4.6 follow from the arguments in Theorems 8 and 14 in [24].

Using simple composition rules and finite induction, it follows from the definition of $\bar{\mathbf{u}}^{K}_{\psi}$ and $\bar{\mathbf {o}}^{K}_{\psi}$ that the continuity of these functions on $I\times P\times\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{x}}$ is a consequence of continuity of u _ψ and o _ψ on $I\times P\times\mathbb{R}^{2n_{x}}\times\mathbb{R}^{6n_{y}}$, provided that $\tilde{\mathbf{y}}$ is chosen as a continuous function on I×P. Similarly, the global Lipschitz condition of Assumption 4.1 (Condition 4) is satisfied provided (u _ψ ,o _ψ )^T is globally Lipschitz with respect to

uniformly on I×P. Since ψ′ is factorable, [24] again provides mild conditions (see Theorem 6) on its factorable form which guarantee continuity of u _ψ and o _ψ on $I\times P\times\mathbb{R}^{2n_{x}}\times \mathbb {R}^{6n_{y}}$, provided that C is continuous on I. If C is Lipschitz on I, then the global Lipschitz condition above follows as well. The definition of ψ′ ensures that the required conditions are satisfied by the factorable representation of ψ′, provided they are satisfied by the factorable representations of g and $\frac{\partial\mathbf{g}}{\partial\mathbf{y}}$ (to assert this, Condition 4 of Assumption 3.2 must be invoked to guarantee that no division by zero occurs in the definition of ψ′). Again, the global Lipschitz condition on u _ψ and o _ψ (Condition 5 of Assumption 4.1) follows from a weaker local condition and the use of the $\operatorname{mid}$ function to map any arguments $(\mathbf {z}_{x}^{c},\mathbf{z}_{x}^{C},\mathbf{z}_{y}^{c},\mathbf{z}_{y}^{C},\tilde{\mathbf {z}}_{y}^{c},\tilde {\mathbf{z}}_{y}^{C},\boldsymbol{\lambda}^{c},\boldsymbol{\lambda}^{C})$ in $\mathbb{R}^{2n_{x}}\times \mathbb{R}^{6n_{y}}$ into the compact interval

in a Lipschitz manner [24].

Appendix B: State Bounds

In this section, the main results in [21, 22] are summarized, culminating in the method used for computing state bounds for (1a)–(1b). The bounds are given as the solution of an auxiliary system of semi-explicit DAEs derived using interval arithmetic and an interval Newton method. To state this system, some preliminaries are required. First, the interval Hansen–Sengupta method is described, which is used for bounding the solutions of nonlinear algebraic equations.

2.1 B.1 The Interval Hansen–Sengupta Method Applied to DAEs

Let $(I,P,Z_{x},Z_{y})\in\mathbb{I}D_{t}\times\mathbb{I}D_{p}\times\mathbb {I}D_{x}\times\mathbb{I}D_{y}$. We are concerned with (i) determining if there exist points z _y ∈Z _y such that g(t,p,z _x ,z _y )=0 for some (t,p,z _x )∈I×P×Z _x , and (ii) computing a refined interval $Z_{y}'\subset Z_{y}$ which contains all such z _y . Conceptually, this is done by using the mean-value theorem. For any $\mathbf{C}\in\mathbb {R}^{n_{y}\times n_{y}}$ and $\tilde{\mathbf{z}}_{y}\in Z_{y}$, consider the interval linear system

(18)

The solution set of (18) is the set of $\boldsymbol{\rho}\in\mathbb{R}^{n_{y}}$ such that A ρ=b for some $\mathbf{A}\in\mathbf {C} [\frac{\partial \mathbf{g}}{\partial\mathbf{y}} ](I,P,Z_{x},Z_{y})$ and $\mathbf {b}\in-\mathbf {C} [\mathbf{g} ](I,P,Z_{x},\tilde{\mathbf{z}}_{y})$. Applying the mean-value theorem as in the proof of Theorem 4.1, one can easily show the following: Any z _y ∈Z _y satisfying g(t,p,z _x ,z _y )=0 for some (t,p,z _x )∈I×P×Z _x must correspond to an element $(\mathbf{z}_{y}-\tilde{\mathbf{z}}_{y})=\boldsymbol{\rho}$ of the solution set of (18). Thus, we are interested in computing an interval enclosure of this solution set.

For Q⊂ℝ, let hull(Q) denote the interval hull of Q; i.e, the smallest interval containing Q. To state the Hansen–Sengupta method formally, the following definition is useful.

Definition B.1

For all $A,B,Z\in\mathbb{IR}$, let

The following lemma provides a way to evaluate Γ computationally. See Proposition 4.3.1 in [25] for proof.

Lemma B.1

For all $A,B,Z\in\mathbb{IR}$,

(19)

where B/A denotes interval division,

For any $A,B,Z\in\mathbb{IR}$, either $\varGamma(A,B,Z)\in\mathbb {IR}$ or Γ(A,B,Z)=∅. For convenience, the definition of Γ is extended, so that Γ(A,B,Z)=∅ when any of A, B, or Z is empty. Furthermore, we adopt the convention that any arithmetic operation between an element of $\mathbb{IR}$ and ∅ returns ∅, and any Cartesian product involving ∅ is equivalent to ∅. The following definition generalizes Γ for application to n dimensional linear systems.

Definition B.2

For $A\in\mathbb{IR}^{n\times n}$, $B,Z\in\mathbb{IR}^{n}$, define Γ(A,B,Z):=W ₁×⋯×W _n , where

Applying Γ to (18) gives the following variant of the well-known result Theorem 5.1.8 in [25], proven in [21].

Corollary B.1

Let $(I,P,Z_{x},Z_{y})\in\mathbb{I}D_{t}\times\mathbb{I}D_{p}\times\mathbb {I}D_{x}\times\mathbb{I}D_{y}$, $\tilde{\mathbf{z}}_{y}\in Z_{y}$, $\mathbf {C}\in\mathbb{R}^{n_{y}\times n_{y}}$ and define

With $Z_{y}':=\mathcal{H}(I,P,Z_{x},Z_{y},\tilde{\mathbf{z}}_{y},\mathbf {C})$, the following conclusions hold:

• 1. If (t,p,z _x ,z _y )∈I×P×Z _x ×Z _y satisfies g(t,p,z _x ,z _y )=0, then $\mathbf{z}_{y}\in Z_{y}'$.

• 2. If $Z_{y}'=\emptyset$, then $\not\exists(t,\mathbf{p},\mathbf {z}_{x},\mathbf {z}_{y})\in I\times P\times Z_{x}\times Z_{y}$ such that g(t,p,z _x ,z _y )=0.

• 3. If $\tilde{\mathbf{z}}_{y}\in\operatorname{int}(Z_{y})$ and $\emptyset \neq Z_{y}'\subset\operatorname{int}(Z_{y})$, then $\exists\mathbf{H}\in C^{1}(I\times P\times Z_{x},Z_{y}')$ such that, for every point (t,p,z _x ) in I×P×Z _x , z _y =H(t,p,z _x ) is the unique element of Z _y satisfying g(t,p,z _x ,z _y )=0. Moreover, the interval matrix $\mathbf{C} [\frac {\partial\mathbf{g}}{\partial\mathbf{y}} ](I,P,Z_{x},Z_{y})$ does not contain a singular matrix and does not contain zero in any of its diagonal elements.

The following theorem is a key result from [21].

Theorem B.1

Let $(I,P,Z_{x},Z_{y})\in\mathbb{I}D_{t}\times\mathbb{I}D_{p}\times\mathbb {I}D_{x}\times\mathbb{I}D_{y}$, $\tilde{\mathbf{z}}_{y}\in Z_{y}$, $\mathbf {C}\in\mathbb{R}^{n_{y}\times n_{y}}$, and define $\mathcal {H}(I,P,Z_{x},Z_{y},\tilde{\mathbf{z}}_{y},\mathbf{C})$ as in Corollary B.1. Furthermore, let $X_{0}\in\mathbb{IR}^{n_{x}}$ satisfy x ₀(P)⊂X ₀ and denote I=[t ₀,t _f ]. If the inclusions

(20)

(21)

(22)

hold, then there exists a regular solution of (1a)–(1b) on I×P satisfying $(\mathbf{x}(t,\mathbf{p}), \mathbf {y}(t,\mathbf{p}))\in Z_{x}\times\ Z_{y}'$ for all (t,p)∈I×P. Furthermore, for any connected set $\tilde{I}\subset I$ containing t ₀, any connected set $\tilde{P}\subset P$, and any solution (x ^∗,y ^∗) of (1a)–(1b) on $\tilde{I}\times\tilde{P}$, either (x ^∗,y ^∗)=(x,y) on $\tilde{I}\times\tilde{P}$, or y ^∗(t ₀,p)∉Z _y , $\forall\mathbf{p}\in\tilde{P}$.

2.2 B.2 Further Definitions

The bounding method of [21, 22] requires some further definitions related to the interval Hansen–Sengupta method.

Definition B.3

Let $\widetilde{\cap}:\mathbb{IR}^{n}\times\mathbb{IR}^{n}\rightarrow \mathbb{IR}^{n}$ be defined by

$$ \widetilde{\bigcap}\bigl(\bigl[\mathbf{z}^L, \mathbf{z}^U\bigr],\bigl[\hat{\mathbf {z}}^L,\hat{\mathbf {z}}^U\bigr]\bigr):=\bigl[\operatorname{mid}\bigl( \mathbf{z}^L,\mathbf{z}^U,\hat{\mathbf {z}}^L \bigr),\operatorname{mid}\bigl(\mathbf{z}^L,\mathbf{z}^U, \hat{\mathbf{z}}^U\bigr)\bigr]. $$

(23)

Furthermore, define the standard notation $Z\,\widetilde{\cap}\,\hat {Z}:=\widetilde{\bigcap}(Z,\hat{Z})$, $\forall Z,\hat{Z}\in\mathbb{IR}^{n}$.

The following definition describes modified forms of Γ used below.

Definition B.4

Let $\mathcal{D}^{*} := \{(A,B,Z)\in\mathbb{IR}^{n\times n}\times \mathbb{IR}^{n}\times\mathbb{IR}^{n}: 0\not\in A_{ii}, \forall i=1,\ldots,n\}$ and define $\varGamma^{*},\varGamma^{+}:\mathcal {D}^{*}\rightarrow\mathbb{IR}^{n}$ by $\varGamma^{*}(A,B,Z):= W^{*}_{1}\times \cdots\times W^{*}_{n}$ and $\varGamma^{+}(A,B,Z):= W^{+}_{1}\times\cdots \times W^{+}_{n_{y}}$, where

(24)

(25)

The following definition formalizes the notation $\mathcal{H}$ from Corollary B.1, with a slight modification to reflect the fact that, in the method presented below, the reference point $\tilde {\mathbf{z}}_{y}$ is taken as the midpoint of Z _y and does not need to be specified independently. Notation is also introduced for iterative application of $\mathcal{H}$, and extended forms based on Γ ⁺ and Γ ^∗ are defined.

Definition B.5

Define $M_{\varGamma}:\mathbb{I}D_{t}\times\mathbb{I}D_{p}\times\mathbb {I}D_{x}\times\mathbb{I}D_{y}\times\mathbb{IR}^{n_{y}\times n_{y}}\rightarrow\mathbb{IR}^{n_{y}\times n_{y}}\times\mathbb {IR}^{n_{y}}\times\mathbb{IR}^{n_{y}}$ by

and define the set

For every K∈ℕ, define the function $\mathcal {H}^{K}:\mathbb{I}D_{t}\times\mathbb{I}D_{p}\times\mathbb{I}D_{x}\times \mathbb{I}D_{y}\times\mathbb{IR}^{n_{y}\times n_{y}}\rightarrow\mathbb {IR}^{n_{y}}$ by $\mathcal{H}^{K}(I,P,Z_{x},Z_{y}^{0},\mathbf{C}):= Z_{y}^{K}$, where $Z_{y}^{k+1}=m(Z^{k}_{y})+\varGamma (M_{\varGamma}(I,P, Z_{x},Z^{k}_{y},\mathbf {C}) )$, for every k in {0,…,K−1}. Furthermore, define $\mathcal{H}^{+,K}:\mathcal{D}_{\mathcal{H}}^{*}\rightarrow \mathbb{IR}^{n_{y}}$ exactly as $\mathcal{H}^{K}$ with Γ ⁺ in place of Γ, and define $\mathcal{H}^{*}:\mathcal{D}_{\mathcal {H}}^{*}\rightarrow\mathbb{IR}^{n_{y}}$ exactly as $\mathcal{H}^{1}$ with Γ ^∗ in place of Γ.

2.3 B.3 State Bounds as the Solutions of an Auxiliary System of DAEs

Definition B.6

Let $\Box:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{IR}^{n}$ be defined by

(26)

Definition B.7

Let $\mathcal{B}_{i}^{L},\mathcal{B}_{i}^{U}:\mathbb{IR}^{n_{x}}\rightarrow \mathbb{IR}^{n_{x}}$ be defined by $\mathcal{B}_{i}^{L}([\mathbf {v},\mathbf {w}])=\{\mathbf{z}\in[\mathbf{v},\mathbf{w}]: z_{i}=v_{i}\}$ and $\mathcal {B}_{i}^{U}([\mathbf{v},\mathbf{w}])=\{\mathbf{z}\in[\mathbf {v},\mathbf{w}]: z_{i}=w_{i}\}$, for every i=1,…,n _x .

We are now prepared to state an auxiliary system of DAEs describing state bounds. Let I=[t ₀,t _f ]⊂D _t , P⊂D _p and X ₀⊂D _x be intervals and suppose that x ₀(P)⊂X ₀. For every i∈{1,…,n _x }, let

(27)

(28)

(29)

(30)

(31)

where

(32)

(33)

(34)

(35)

(36)

Choosing any K∈ℕ, define the functions in (27)–(31) by

(37)

(38)

(39)

(40)

(41)

(42)

(43)

(44)

(45)

For any continuous and pointwise positive γ:I→ℝ, consider the initial value problem in DAEs

(46)

(47)

(48)

(49)

for all i=1,…,n _x , with initial conditions

(50)

In (48) and (49), $\mathcal{H}^{*,L}$ and $\mathcal{H}^{*,U}$, denote the lower and upper bounds of $\mathcal{H}^{*}$, respectively. A solution of (46)–(50) is any $\mathbf {v},\mathbf{w}:I\rightarrow\mathbb{R}^{n_{x}}$ and $\mathbf {z}_{y}^{L},\mathbf {z}_{y}^{U}:I\rightarrow\mathbb{R}^{n_{y}}$ satisfying (46)–(50) for all t∈I, with v and w continuously differentiable and $\mathbf {z}_{y}^{L}$ and $\mathbf{z}_{y}^{U}$ piecewise C ¹ (see [35] for a definition of this class of functions). The main result of [22], stated below, shows that any such solution provides state bounds for a unique regular solution of (1a)–(1b) on I×P. The results below require that the interval extensions [f], [g] and $[\frac{\partial\mathbf {g}}{\partial\mathbf {y}}]$ are piecewise C ¹, which is defined and discussed in detail in [22]. In the remainder of this section, we use the notation $Z_{y}(t):=[\mathbf{z}_{y}^{L}(t),\mathbf{z}_{y}^{U}(t)]$ and $Z_{y}'(t):=\mathcal {H}^{*}(\phi(t,\mathbf{v}(t),\mathbf{w}(t),\mathbf {z}_{y}^{L}(t),\mathbf{z}_{y}^{U}(t)))$.

Lemma B.2

Let $(\mathbf{v},\mathbf{w},\mathbf{z}_{y}^{L},\mathbf{z}_{y}^{U})$ be a solution of (46)–(50) on I and define

Then

• 1. $[\mathbf{v}(t),\mathbf{w}(t)]\times Z_{y}(t)\in\mathbb {I}D_{x}\times \mathbb{I}D_{y}$, ∀t∈I.

• 2. There exists H∈C ¹(V,D _y ) such that, for every (t,p,z _x )∈V, z _y =H(t,p,z _x ) is an element of $Z_{y}'(t)$ and satisfies g(t,p,z _x ,z _y )=0 uniquely among elements of Z _y (t).

• 3. For every t∈I, the interval matrix

  

  does not contain any singular matrix and does not contain zero in any of its diagonal elements.

Theorem B.2

Let $(\mathbf{v},\mathbf{w},\mathbf{z}_{y}^{L},\mathbf{z}_{y}^{U})$ be a solution of (46)–(50) on I. Then there exists a regular solution (x,y) of (1a)–(1b) on I×P satisfying

(51)

(52)

for all (t,p)∈I×P and any q∈ℕ. Furthermore, for any connected $\tilde{I}\subset I$ containing t ₀, any connected $\tilde{P}\subset P$, and any solution (x ^∗,y ^∗) of (1a)–(1b) on $\tilde{I}\times\tilde{P}$, either (x ^∗,y ^∗)=(x,y) on $\tilde {I}\times\tilde {P}$, or y ^∗(t ₀,p)∉Z _y (t ₀), $\forall \mathbf{p}\in \tilde{P}$.

In light of Theorem B.2, state bounds are given by simply solving the DAEs (46)–(50). Provided that numerical error is not a critical concern, this can be done using any state-of-the-art DAE solver. For the numerical example in Sect. 5, the code IDA [36] was used with absolute and relative tolerances of 10⁻⁶. Furthermore, K=5 and γ(t)=10⁻⁵, ∀t∈I.