Lars Grüne

Manuela Sigurani

Numerical ISS controller design via a dynamic game approach

We present a numerical design method for an input-to-state practically stabilizing (ISpS) feedback controller for perturbed discrete time nonlinear control systems. By appropriately scaling the system, the problem of ISpS controller design can be converted into a uniform practical stabilization problem under perturbations which can be solved by a numerical dynamic game approach. We present and analyze this approach with a particular focus on a quantitative analysis of the resulting gain and the size of the exceptional region for practical stability.

I. INTRODUCTION

Since its introduction by Sontag in [START_REF] Sontag | Smooth stabilization implies coprime factorization[END_REF], the concept of input-to-state stability has become one of the most influential concepts in nonlinear stability under perturbations. Besides yielding a theoretically sound concept for the quantitative and qualitative analysis of stability of nonlinear systems under perturbations, one of its particular features is the possibility to analyze the stability of interconnected systems by means of analyzing low dimensional subsystems via ISS small gain arguments, see, e.g., [START_REF] Jiang | Small-gain theorem for ISS systems and applications[END_REF], [START_REF] Dashkovskiy | An ISS small gain theorem for general networks[END_REF], [START_REF] Dashkovskiy | Small gain theorems for large scale systems and construction of ISS Lyapunov functions[END_REF]. The latter particularly allows for a rigorous approach to decentralized controller design by designing input-to-state stabilizing controllers. Hence, a systematic way to design ISS controllers can in particular serve as a building block for decentralized stabilizing feedback design for large networks of nonlinear systems.

Similar to stabilizing controllers for nonlinear systems, ISS controllers can in principle be derived from corresponding ISS Lyapunov functions via a universal or Sontag type formula [START_REF] Liberzon | Universal construction of feedback laws achieving ISS and integral-ISS disturbance attenuation[END_REF]. However, applying this formula requires the analytic knowledge of an ISS Lyapunov function which may not always be available. An alternative are dynamic programming type design methods relying on optimal control formulations whichrequire much less a priory analytical knowledge. The main drawback of this approach is the curse of dimensionality which means that it is only computationally feasible for systems of moderate space dimension. However, in a small gain based decentralized setting in which the controller design is to be carried out for a set of low dimensional subsystems the approach may provide an attractive and feasible alternative, which is why we investigate it in this paper.

A dynamic programming based ISS controller design was proposed in [START_REF] Huang | A unified approach to controller design for achieving ISS and related properties[END_REF] 1 . However, the approach in this paper has the drawback that by converting the problem into an auxiliary ℓ ∞ control problem the state variable needs to be augmented by 2 additional scalar states which considerably increases the computational complexity of the controller design. In order to avoid this problem, in this paper we propose an approach which consists of converting the ISS controller design problem into a uniform stabilization problem for a perturbed system which can be accomplished without increasing the dimension. The equivalence between ISS and robust stability was already exploited in a theoretical context in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF] (for a continuous time version of this result see [START_REF] Sontag | New characterizations of input-tostate stability[END_REF]) and thus our approach can be seen as a constructive numerical interpretation of the results in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF]. In order to solve the auxiliary stabilization problem under perturbations we use the game theoretic algorithmic approach from [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] which in turn relies on [START_REF] Junge | Global optimal control of perturbed systems[END_REF], [START_REF] Osinga | A set oriented approach to global optimal control[END_REF]. Like in most nonlinear numerical approaches relying on Lyapunov functions, cf. e.g., [START_REF] Giesl | Construction of global Lyapunov functions using radial basis functions[END_REF], [START_REF] Hafstein | An algorithm for constructing Lyapunov functions[END_REF], [START_REF] Giesl | Existence of piecewise linear Lyapunov functions in arbitrary dimensions[END_REF], a neighborhood of the equilibrium (in this paper always chosen as the origin) needs to be treated in a different way. In our setting this means that the resulting nonlinear controller will in general only yield input-to-state practical stability (ISpS), i.e., the closed loop system will have the ISS property only outside a neighborhood of the origin. For this reason, a substantial part of the analysis in this paper is devoted to keeping track of the size of this exceptional neighborhood in order to control the errors introduced by the numerical solution of the dynamic game problem. This allows us to identify conditions under which this neighborhood is small. Hence, a linearization based design could be used if desired in order to define an ISS controller also near the origin. Moreover, we note that both the Lyapunov function as well as the resulting optimal feedback law are piecewise constant and thus discontinuous in our approach, which is why we provide an analysis entirely avoiding continuity assumptions.

II. SETTING

Our goal is to construct an input-to-state practically stabilizing (ISpS) controller for the controlled and perturbed discrete-time system

x k+1 = f (x k , u k , w k ), k = 0, 1, . . . (1)
For simplicity of exposition, in this paper we consider a compact state space X ⊂ R d . Extensions to hybrid systems with state spaces containing discrete state components and more general set stabilization problems can be addressed similarly to [START_REF] Junge | Optimal stabilization of hybrid systems using a set oriented approach[END_REF]. The discrete time model under consideration can, of course, be the discrete time representation of a sampled continuous time model. The values u k and w k denote the control and perturbation acting on the system which are taken from sets U ⊂ R m an W ⊂ R q , respectively, which again for simplicity of exposition are supposed to be compact. Infinite sequences of control and perturbation values are denoted by u = (u 0 , u 1 , . . .) and w = (w 0 , w 1 , . . .) and the corresponding spaces of such sequences with values u k ∈ U and w k ∈ W are denoted by U and W , respectively.

The control objective of designing an ISpS controller means that we intend to find a static state feedback controller u k = u P (x k ) such that the closed loop system

x k+1 = f (x k , u P (x k ), w k ) =: h(x k , w k ), k = 0, 1, . . . (2)
is input-to-state practically stable in the following sense.

Definition 1: System ( 2) is called input-to-state practically stable (ISpS) with respect to δ , ∆ w ∈ R ≥0 on a set Y ⊂ X if there exist β ∈ K L and γ ∈ K , such that the solutions of the system satisfy

x k ≤ max { β ( x 0 , k), γ( w ∞ ), δ } , (3)
for all x 0 ∈ Y , all w ∈ W with w ∞ ≤ ∆ w and all k ∈ N 0 .

The approach we present in this paper relies on the conversion of the ISpS controller design method into a uniformly practically stabilizing controller design method. To this end, in the next section we first sketch the dynamic game based uniformly practically stabilizing controller design method from [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF]. Afterwards, we explain how to use this approach for the ISpS controller design problem at hand.

III. GAME THEORETIC STABILIZING CONTROLLER DESIGN FOR PERTURBED SYSTEMS

In this section we consider the perturbed control system

x k+1 = f (x k , u k , d k ), k = 0, 1, . . . (4)
Except for changing the notation, this model has the same structure as [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]. Moreover, x k , u k and the respective sets and spaces remain unchanged. The perturbation values d k are now taken from a set D ⊂ R q , the corresponding sequences are denoted as d = (d 0 , d 1 , . . .) and the space of such sequences with d k ∈ D is denoted by D.

The control objective for System (4) is to design a practically uniformly stabilizing state feedback controller, i.e., a controller u k = u P (x k ) such that the closed loop system

x k+1 = f (x k , u P (x k ), d k ) =: h(x k , d k ), k = 0, 1, . . . (5)
satisfies the following definition.

Definition 2: System (5) is called uniformly (w.r.t. d ∈ D) practically (w.r.t. δ ) asymptotically stable on a set Y ⊂ X if there exist β ∈ K L such that the solutions of the system satisfy

x k ≤ max { β ( x 0 , k), δ } , (6)
for all x 0 ∈ Y , all d ∈ D and all k ∈ N 0 .

In order to design a controller we employ the dynamic game approach from [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] which in turn relies on ideas from [START_REF] Junge | Global optimal control of perturbed systems[END_REF], [START_REF] Osinga | A set oriented approach to global optimal control[END_REF]. For a given initial state x ∈ X, a given control sequence u = (u k ) k∈N ∈ U and a given perturbation sequence d = (d k ) k∈N ∈ D, we denote the solution trajectory of ( 4) by

(x k (x, u, d)) k∈N .
Specifying a target set T (typically a small neighborhood of the origin) and a stage cost g(x, u), we define the accumulated cost as

J(x, u, w) = k(T,x,u,d) ∑ k=0 g(x k (x, u, d), u k ), with k(T, x, u, d) := inf{k ≥ 0 | x k (x, u, d) ∈ T }.
Here we assume that the stage cost satisfies the following assumption.

Assumption 3: The stage cost g penalizes the distance to 0, i.e., there exists α ∈ K ∞ such that

g(x, u) ≥ α( x ) (7)
holds for all x ∈ X, u ∈ U.

Using the accumulated cost we now define a dynamic game by means of its upper value function

V (x) = sup β ∈B inf u∈U J(x, u, β (u)), x ∈ X, (8)
Here, B denotes the set of all nonanticipating strategies β :

U → D, i.e. all strategies β satisfying u k = u ′ k ∀k ≤ K ⇒ β (u) k = β (u ′ ) k ∀k ≤ K for any two control sequences u = (u k ) k , u ′ = (u ′ k ) k ∈ U .
By standard dynamic programming arguments [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF] one sees that this function fulfills and is uniquely determined by the optimality principle

V (x) = inf u∈U sup d∈D g(x, u) +V ( f (x, u, d)) (9)
for x / ∈ T together with the boundary condition V | T ≡ 0. This equation can also be written as

V (x) = inf u∈U g(x, u) + sup x ′ ∈ f (x,u,D) V (x ′ ) . (10)
Note that in [START_REF] Üller | Set oriented optimal control using past information[END_REF] the parameterization of f by d is not needed any more since it is sufficient to know the set valued image f (x, u, D). This fact is exploited in the discretization proposed in [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF]. In this discretization, the set X is decomposed into a finite partition P of boxes or cells P with pairwise disjoint interior and P∈P P = X. Then, a graph theoretic representation of the dynamics on P is constructed. Since the model includes both control and perturbation, the resulting graph theoretic approximation takes the form of a hypergraph. Solving a generalized min-max shortest path problem on this hypergraph then yields an approximation V P of V which is constant on each element P of the partition P.

The central trick introduced in [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] in terms of stabilization is to interpret the discretization error introduced by the partition P as a perturbation and explicitly include it in the computation. This way, one can prove that the approximation V P is not an arbitrary approximation to V but that V P is in fact a practical control Lyapunov function. While in [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] the discretization error is the only perturbation acting on the system, here we extend the setting by considering both the original perturbation d and the discretization error as perturbations. The necessary changes in the statements and proofs of the results in [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] are straightforward and are hence omitted due to space constraints. The resulting theorem then reads as follows.

Theorem 4: ([8, Theorem 1]) Let V denote the optimal value function of the dynamic game (8) and let V P denote the approximate optimal value function constructed according to the algorithm presented in [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] on a given partition P with target set T ⊂ P and 0 ∈ T . Then,

V (x) -max y∈T V (y) ≤ V P (x), (11)
i.e. V P is an upper bound for V -maxV | T . Furthermore, V P satisfies

V P (x) ≥ min u∈U g(x, u) + sup x ′ ∈ f (x,u,D) V P (x ′ ) (12)
for all x ∈ X \ T . Note that V P may assume the value +∞ on some parts of X, in which case inequality [START_REF] Ne Ši Ć | Optimization based stabilization of sampled-data nonlinear systems via their approximate discrete-time models[END_REF] does not yield valuable information. We define the stabilizable set w.r.t. V P by

S P := {x ∈ X |V P (x) < ∞}.
Theorem 4 now motivates the definition of the controller

u P (x) := argmin u∈U g(x, u) + sup x ′ ∈ f (x,u,D) V P (x ′ ) (13)
for x ∈ S P \ T . For x ∈ T we may either use Formula (13) or we may define u P in a different way, e.g., constantly equal to 0 or by using linearization at the origin, as, e.g., proposed in [START_REF] Üne | Two complementary approaches to event-based control, at-Automatisierungstechnik[END_REF] in an event based context. For x ∈ X \S P , our approach does not allow for a meaningful definition of u P .

The following theorem from [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] summarizes the properties of u P .

Theorem 5: ([8, Theorem 2]) Consider system (4), a target set T ⊂ P with 0 ∈ T and the optimal value function V P constructed according to the algorithm presented in [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] on a given partition P. For x 0 ∈ S P denote the trajectory of the closed loop system (5) with feedback u = u P from (13) by x k .

Then there exists k * ∈ N such that

V P (x k ) ≥ g(x k , u P (x k )) +V P (x k+1 ), k ∈ {0, . . . , k * -1}, x k * ∈ T and x k ∈ N, k > k *
where N := {x ∈ X : V P (x) ≤ ν} and ν := max x∈T,d∈D V P ( f (x, u P (x), d)), i.e., the closed loop trajectory enters the target set T at time k * and stays in the sublevel set N, afterwards. Remark 6: Note that it is in general not possible to use the target set T = {0} unless we are willing to assume that the system can be controlled to the origin on finitely many steps (and even then using T = {0} is likely to cause numerical problems). Similar problems in small neighborhoods of the equilibrium occur in many other numerical approaches for computing Lyapunov functions for nonlinear systems, even for non-controlled systems, see [START_REF] Giesl | Construction of global Lyapunov functions using radial basis functions[END_REF], [START_REF] Hafstein | An algorithm for constructing Lyapunov functions[END_REF], [START_REF] Giesl | Existence of piecewise linear Lyapunov functions in arbitrary dimensions[END_REF]. This means that on a small neighborhood around the origin V P is not a Lyapunov function and consequently the exceptional set N cannot be avoided. The size of this set then determines the size of the parameter δ on the practical stability definition.

In control problems, the usual way to work around this problem is to use linearization techniques in order to solve the feedback stabilization problem locally near 0. For this purpose it is of utmost importance to keep the size of the exceptional set N and thus of δ small. Consequently, one of the central tasks in the following section will be to carefully estimate this value in the ISpS context.

IV. ISPS CONTROLLER DESIGN

In order to apply the algorithm from the previous section to practical ISS controller design we make use of one of the central results in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF], which states that System (2) is ISS if and only if it is robustly stable, i.e., if there exists e : R n × R q → R q and η ∈ K ∞ such that System (5) with

f (x, u, d) = f (x, u, e(x, d)) and D = B 1 (0) (14)
is uniformly asymptotically stable, where e is such that for each w ∈ W with w ≤ η( x ) there exists d ∈ D with e(x, d) = w. For instance, e could be defined as e(x, d) := η( x )d which is also the choice in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF]. The proof of the equivalence between ISS and robust stability relies on Lyapunov function arguments.

In order to carry over the proof to our setting of practical stability, in this section we first provide a characterization of ISpS by means of an appropriate Lyapunov function V . In contrast to the analogous non-practical statements in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF], here we give a direct proof which allows to determine the resulting gains and the size of the practical stability region. Afterwards we then show that V = V P is an ISpS Lyapunov function in our sense for the non-scaled system.

A particular difficulty in this derivation is the fact that in general neither V P nor u P are continuous. Hence, we will need to avoid all arguments in the proof which require continuity of V P or u P . We begin by introducing a continuity condition on f and a boundedness condition on h.

Assumption 7: The map f : 1) is uniformly continuous in w in the following sense: there exist

X × U × W → R n in (
γ w ∈ K ∞ , such that for all x ∈ X, u ∈ U and w ∈ W f (x, u, w) -f (x, u, 0) ≤ γ w ( w ). (15)
Since we assume X and U to be compact it is not very restrictive to assume that γ w is independent of x and u. If needed, the inequality could be suitably generalized.

Assumption 8: The map h : X ×W → R n in (2) is bounded for w = 0 on T in the following sense: there exists γ x ∈ K ∞ such that for all sufficiently small sets T ⊂ X with 0 ∈ T and each x ∈ T we have

h(x, 0) ≤ γ x ( x ). (16)
Observe that Assumption 8 is satisfied, e.g., if u P ≡ 0 on T and f (x, 0, 0) is Lipschitz on T with constant L. In this case, since h(0, 0) = f (0, 0, 0) = 0, we may choose γ x (r) = Lr. Definition 9: A continuous function V : X → R ≥0 is called an ISpS Lyapunov function for system (2) on a sublevel set Y = {x ∈ X |V (x) ≤ ℓ} for some ℓ > 0 if there exist functions α, α ∈ K ∞ , α, ρ ∈ K , a value w ∈ R >0 , and a set T ⊂ Y such that for all x ∈ Y \ T the inequalities and implications

α( x ) ≤ V (x) ≤ α( x ) (17)
and

V (x) ≥ ρ( w ) ⇒ V (h(x, w)) -V (x) ≤ -α( x ) (18)
hold for all w ∈ W with w ≤ w.

The next theorem now shows that the existence of an ISpS Lyapunov function implies ISpS and specifies the resulting gain γ and the size δ in the practical stability definition.

Theorem 10: Consider system (2) satisfying Assumptions 7 and 8 and assume that the system admits an ISpS Lyapunov function V . Let a = sup x∈T x . Then the system is ISpS on Y with

δ = α -1 • α max{2γ x (a), a} , γ(r) = α -1 • α max{2γ w (r), α -1 (ρ(r)) + γ w (r)} and ∆ w = γ -1 (α -1 (ℓ)) provided δ ≤ α -1 (ℓ) holds.
Proof: We fix x 0 ∈ Y , w ∈ W and denote the corresponding trajectory of (2) by x k . We begin the proof by deriving estimates for V (x k ) under different assumptions. To this end, we distinguish four different cases.

Case 1: Let k ′ ∈ N be such that V (x k ) ≥ ρ( w ∞ ) and x k ∈ Y \ T for all k = 0, ..., k ′ -1. Then (17) yields

x(k) ≥ α -1 (V (x(k))) (19)
and from ( 18) we obtain

V (x k+1 ) -V (x k ) (18) ≤ -α( x k ) ( 20
) (19) ≤ -α α -1 (V (x k )) . (21)
Note that x 0 ∈ Y and the definition of Y implies x k ∈ Y for all k = 0, . . . , k ′ -1, hence (18) may indeed be used for all these k. Setting α := α • α -1 , Lemma 4.3 from [START_REF] Jiang | A converse Lyapunov theorem for discrete-time systems with disturbances[END_REF] then yields the existence of β ∈ K L such that

V (x(k)) ≤ β (V (x 0 ), k) for all k = 0, . . . , k ′ . ( 22
)
Case 2:

Let k ∈ N be such that V (x k ) < ρ( w ∞ ) and x k ∈ Y \ T . This implies α( x k ) ≤ V (x(k)) < ρ( w ∞ ). ( 23
)
and since [START_REF] Jiang | A converse Lyapunov theorem for discrete-time systems with disturbances[END_REF] implies

V (h(x k , 0)) ≤ V (x k ) we get h(x k , 0) ≤ α -1 (ρ( w ∞ ).
Then Assumption 7 implies

x k+1 ≤ h(x k , 0) + γ w ( w ∞ ) ≤ α -1 (ρ( w )) + γ w ( w ∞ ) =: γ( w ∞ ).
This yields

V (x k+1 ) ≤ α(x k+1 ) ≤ α( γ( w ∞ )) =: γ( w ∞ ). Case 3: Let k ∈ N be such that V (x k ) < γ( w ∞ ) and x k ∈ Y \ T . Then we either have V (x k ) < ρ( w ∞ ) and thus Case 2 implies V (x k+1 ) ≤ γ( w ∞ ).
Otherwise, we have V (x k ) ≥ ρ( w ) and ( 18) yields

V (x k+1 ) ≤ V (x k ) ≤ γ( w ∞ ).
Thus, in either case we get V (x k+1 ) ≤ γ( w ∞ ). Case 4: Let k ∈ N be such that x k ∈ T . Then, from Assumptions 7 and 8 we get

x k+1 = h(x k , w k ) ≤ h(x k , 0) + h(x k , w k ) -h(x k , 0) ≤ γ x ( x k ) + γ w ( w ∞ ).
Defining a := sup x∈T x we can continue

x k+1 ≤ γ x (a) + γ w ( w ∞ ) ≤ max{2γ x (a), 2γ w ( w ∞ )}
and using ( 17) we get

V (x k+1 ) ≤ α(max{2γ x (a), 2γ w ( w ∞ )}) = max{α(2γ x (a)), α(2γ w ( w ∞ ))}.
Combining these four cases we can now prove the desired inequality (3):

Let k ′ ∈ N be maximal such that the condition from Case 1 is satisfied. Then, for all k = 0, . . . , k ′ we get

x k (17) ≤ α -1 (V (x k )) (22) ≤ α -1 ( β α (V (x 0 ), k)) (17) ≤ α -1 ( β α ( α( x 0 ), k)) =: β ( x 0 , k) (24)
Now, for all k ≥ k ′ by induction we show the inequality

V (x k ) ≤ max{α(2γ x (a)), α(a), α(2γ w ( w ∞ )), γ( w ∞ )}.
(25) Note that the bounds on δ and ∆ w in the assertion ensure that (25) implies V (x k ) ≤ ℓ and thus x k ∈ Y for all w ∈ W with w ∞ ≤ ∆ w . Hence, (25) implies that one of the Cases 1-4 must hold for x k . Consequently, if we know that (25) holds we can use the estimates in the Cases 1-4 in order to conclude an inequality for x k+1 .

To start the induction at k = k ′ , one checks that by definition γ(r) ≥ ρ(r) holds, hence (25) follows from the condition in Case 1 and the maximality of k ′ .

For the induction step k → k + 1, assume that (25) holds for x k . Then, either Case 1 holds implying V (x k+1 ) ≤ V (x k ) and thus (25) for x k+1 . Otherwise, one of the Cases 2-4 must hold for x k which also implies (25) for x k+1 .

Together, (24) and (25) show that either x k ≤ β ( x 0 , k) or x k ≤ max{γ( w ∞ ), δ } holds, which shows the desired ISpS inequality [START_REF] Dashkovskiy | Small gain theorems for large scale systems and construction of ISS Lyapunov functions[END_REF].

The following proposition shows that V P when computed from (4), ( 14) is an ISpS Lyapunov function for [START_REF] Dashkovskiy | An ISS small gain theorem for general networks[END_REF].

Proposition 11: Consider the system (1), the function V P satisfying Theorem 4 for system (4) with f from ( 14) and the corresponding feedback u P from [START_REF] Pannek | Nonlinear Model Predictive Control. Theory and Algorithms[END_REF]. Then V P is an ISpS Lyapunov function for the closed loop system (2) for any ℓ > 0.

Proof: Proof of (17): Obviously, V P (x) ≥ inf u∈U g(x, u) if x / ∈ T . Hence the lower bound in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF] follows directly from Assumption 3. The existence of an upper bound follows since V P ≡ 0 holds on T and V P is piecewise constant and bounded by ℓ on Y . Hence, sup x∈Y, x ≤r V P (x) is piecewise constant, finite for each r > 0 and equal to 0 for all sufficently small r > 0. Thus, it can be overbounded by a function α ∈ K ∞ (which could, e.g., be constructed by piecewise linear interpolation).

Proof of (18): From ( 12) and the definition of u P via (13) it follows that V P (x) ≥ g(x, u P (x)) +V P ( h(x, d)) holds for all x ∈ Y \ T and all d ∈ D. Since by assumption on e in ( 14) for all w ∈ W with w ≤ η( x ) we find d ∈ D with w = e(x, d), this implies that for all such w we obtain the right inequality in [START_REF] Jiang | A converse Lyapunov theorem for discrete-time systems with disturbances[END_REF]. This shows that [START_REF] Jiang | A converse Lyapunov theorem for discrete-time systems with disturbances[END_REF] holds with ρ = η.

Note that since V P assumes only finitely many different values and is finite on S P , choosing ℓ := max x∈S P V P (x) yields the maximal possible domain Y = S P on which V P is an ISpS Lyapunov function.

Our final and main result now summarizes the conditions under which the feedback u P indeed renders System (1) ISpS.

Theorem 12: Consider a system (1) satisfying Assumptions 7 and 8, the function V P satisfying Theorem 4 for system (4) with f from ( 14) and the corresponding feedback u P from [START_REF] Pannek | Nonlinear Model Predictive Control. Theory and Algorithms[END_REF]. Let ℓ ≤ max x∈S P V P (x) and let α, α ∈ K ∞ denote functions such that (17) holds2 on Y = {x ∈ X |V P (x) ≤ ℓ}.

(i) If the value ℓ > 0 is such that the inequality

ℓ ≥ α max{2γ x (a), a} (26)
for a = sup x∈T x holds, then the system is ISpS on Y w.r.t δ = α -1 • α max{2γ x (a), a} and ∆ w as specified in Theorem 10.

(ii) If for each target T which forms a neighborhood of 0 there exists a partition P such that α in (17) can be chosen independently of T , then for each δ > 0 there exists T and P such that u P is ISpS on Y w.r.t. this δ and ∆ w as specified in Theorem 10.

Proof: (i) By Proposition 11 the function V P is an ISpS Lyapunov function. Since (26) ensures that Theorem 10 is applicable, this yields the ISpS property.

(ii) By choosing T to be a sufficiently small neighborhood around 0 we can choose a = sup x∈T x arbitrarily close to 0. Since, moreover, α is independent of the choice of T (and by the first part of the proof of Proposition 11 α is independent of T , too) we can always ensure that (26) holds and that δ as specified in (i) can be chosen arbitrarily small. This shows the assertion.

Remark 13: Under suitable conditions it is reasonable to expect that the condition on α from Theorem 12(ii) holds provided the partition P is -at least near the target Tsufficiently fine. This is due to the fact that with similar arguments as in the proofs of [START_REF] Pannek | Nonlinear Model Predictive Control. Theory and Algorithms[END_REF]Theorem 4.3] or [START_REF] Ne Ši Ć | Optimization based stabilization of sampled-data nonlinear systems via their approximate discrete-time models[END_REF]Theorem 5.4] we can first ensure that such a bound holds for V . In the next step, using similar arguments as in [START_REF] Junge | Global optimal control of perturbed systems[END_REF], by choosing the partition P sufficiently fine one can ensure that V P is close to V from which the assertion follows. Note that the smaller the target T is the smaller the error and thus the finer the partition P must be in order to make this argument work. Details about this construction will be investigated in future research.

Remark 14: Since the computational part of our approach entirely relies on computing a uniformly practically asymptotically stabilizing feedback law for the scaled system (4) by means of (a slight generalization of) the approach from [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF], all the extensions developed for this algorithm can be readily applied to our proposed approach. Particularly, this concerns the approach to implement the computation of feedback laws depending not only on the current but also on past values of the state [START_REF] Üller | Set oriented optimal control using past information[END_REF] and an event based computation of u P [START_REF] Üller | An algorithm for event-based optimal feedback control[END_REF]. Both extensions together allow for a significant reduction of the number of partition elements representing u P . They have thus been used in the computations of our numerical example in the next section.

V. NUMERICAL EXAMPLE

In order to illustrate our approach we show numerical results for a thermofluid process. The process consists of the fill level x 1 and the temperature x 2 of a liquid in a tank. The inflow of liquid can be controlled by u 1 and the liquid in the tank can be cooled using u 2 . The perturbations of the water level w 1 and the temperature w 2 model the unknown inflow of liquid from a second tank. For a description of the full model see [START_REF] Lunze | Control Theory of Digitally Networked Systems[END_REF]Appendix]. After some simplifications of the equations, the behavior of the tank system is described by the state-space model

ẋ1 (t) = 1 0.065 161 • 10 -6 u 1 (t) + 129 • 10 -6 w 1 (t) + 0.34 -270 • 10 -6 x 1 (t) ẋ2 (t) = 1 0.065x 1 (t) 129 • 10 -6 w 1 (t) + 0.34 × (w 2 (t) + 300 -x 2 (t)) +97 • 10 -6 u 2 (t) (287 -x 2 (t)) with X = [ 0.25, 0.4 ] × [ 290, 320 ], w 1 ∈ [ -0.09, 0.09 ], w 2 ∈ [ -20, 20 ] and u i ∈ [ 0, 1 ], i = 1, 2.
For u ⋆ = (0.481465, 0.48466) T , the equation exhibits the equilibrium x ⋆ = (0.32, 295) T . Note that x ⋆ is asymptotically stable, hence the goal of our ISpS controller is not to stabilize the system at x ⋆ but to increase the robustness of the stability against perturbations.

As the system has a cascaded (or triangular) structurei.e., the first equation does not depend on x 2 -it has turned out beneficial to choose e in [START_REF] Hafstein | An algorithm for constructing Lyapunov functions[END_REF] to reflect this structure, i.e., to have the first component independent of x 2 . Note that this way the assumption needed in the second part of the proof of Proposition 11 is not satisfied; however, the proof can be adapted to the cascaded situation. Hence we chose e in ( 14) as

e(x, d) = 1.25297(x 1 -x ⋆ 1 ) 2 d 1 618.75(x 1 -x ⋆ 1 ) 2 + 0.6273(x 2 -x ⋆ 2 ) 2 d 2 We computed the controller using the stage cost g(x, u) = 4 • 10 4 (x 1 -x ⋆ 1 ) 2 + (x 2 -x ⋆ 2 )
2 on a partition P consisting of 32×32 equally sized elements where the target T was chosen as the partition element containing x ⋆ . The sampling time was chosen as 2s, and the control and perturbation value sets U = [0, 1] 2 and D = [-1, 1] 2 were discretized with grids of 9 × 5 and 3 × 3 equidistant nodes, respectively.

For the trajectory simulations a randomly generated sequence w of perturbations was used, using uniformly distributed random numbers in [0, 0.09] and [0, 20], respectively, for the components of each vector w k ∈ R 2 . The resulting trajectories with and without control (both for the same sequence w) are shown in Figure 1. One clearly sees that the controller is able to bring the system considerably closer to the desired equilibrium. However, the practical nature of the controller is also clearly visible by the zig-zagging effect of the x 1 -component. This effect could probably be reduced by using a local robust controller near x ⋆ instead of the constant equilibrium control value u ⋆ we have employed in our simulation.

VI. CONCLUSIONS AND OUTLOOK

We have shown that the equivalence between ISS and robust stability proved in [START_REF] Jiang | Input-to-state stability for discrete-time nonlinear systems[END_REF] for discrete time systems can be used as the basis for a numerical approach for computing ISS feedback laws. As the underlying algorithm for calculating uniformly stabilizing feedback laws from [START_REF] Junge | Approximately optimal nonlinear stabilization with preservation of the Lyapunov function property[END_REF] yields only practically stabilizing controllers, the resulting feedback law will be input-to-state practically stabilizing w.r.t. some δ > 0. In order to obtain a meaningful stability property, a careful analysis of the size δ of the practical stability region is provided in this paper. As Theorem 12 reveals, the existence of an upper bound α on V P which is independent of the target set T is a crucial property for bounding δ .

The rigorous derivation of such upper bounds will be investigated in future research, cf. Remark 13. Moreover, we intend to address the extension of our analysis to event based ISpS feedback laws and the proposed computational approach for ISpS controllers will be used as a building block for a distributed feedback design for large networks of systems based on the small gain arguments from [START_REF] Dashkovskiy | An ISS small gain theorem for general networks[END_REF], [START_REF] Dashkovskiy | Small gain theorems for large scale systems and construction of ISS Lyapunov functions[END_REF].

2 Fig. 1 .

 21 Fig. 1. Trajectories with control (solid) and without control (dashed), x 1component (left) and x 2 -component (right)

This work was supported by the DFG Priority Program 1305 and by the European Union under the 7th Framework Programme FP7-PEOPLE-2010-ITN, Grant agreement nr. 264735-SADCO L. Grüne and M. Sigurani are with the Mathematical Institute, University of Bayreuth, 95440 Bayreuth, Germany lars.gruene, manuela.sigurani@uni-bayreuth.de

This reference treats various different robust controller design objectives, among them ISS.

These functions exist according to the first part of the proof of Proposition 11.


Sigurður Hafstein

Christopher M Kellett

Huijuan Li

Continuous and Piecewise Affine Lyapunov Functions using the Yoshizawa Construction

We present a novel numerical technique for the computation of a Lyapunov function for nonlinear systems with an asymptotically stable equilibrium point. Our proposed approach constructs a continuous piecewise affine (CPA) function given a suitable partition of the state space, called a triangulation, and values at the vertices of the triangulation. The vertex values are obtained from a Lyapunov function in a classical converse Lyapunov theorem and verification that the obtained CPA function is a Lyapunov function is shown to be equivalent to verification of several simple inequalities. Furthermore, by refining the triangulation, we show that it is always possible to construct a CPA Lyapunov function. Numerical examples are presented demonstrating the effectiveness of the proposed method.

I. INTRODUCTION

Lyapunov's Second or Direct Method [START_REF] Lyapunov | The general problem of the stability of motion[END_REF] (see also [START_REF] Rouche | Stability Theory by Liapunov's Direct Method[END_REF]) has proved to be one of the most useful tools for demonstrating stability properties. This is largely due to the fact that if one has a Lyapunov function at hand there is no need to explicitly generate system solutions in order to determine stability. Unfortunately, this frequently trades the difficult problem of generating system solutions for the equally difficult problem of constructing a Lyapunov function. Nonetheless, it is frequently worthwhile numerically searching for a Lyapunov function as it can be used for purposes other than stability analysis such as in estimating a basin of attraction (in the case of local stability) or for inferring robustness properties [START_REF] Kellett | On the robustness of KL-stability for difference inclusions: Smooth discrete-time Lyapunov functions[END_REF] [START_REF] Teel | A smooth Lyapunov function from a class-KL estimate involving two positive semidefinite functions[END_REF].

So-called converse Lyapunov theorems provide existence results for Lyapunov functions; i.e., assuming a particular stability property holds then there exists an appropriate Lyapunov function. However, such results are largely not constructive in nature and, in fact, depend explicitly on solutions of the system under study. As a consequence, various approaches have been proposed for the numerical construction of Lyapunov functions such as collocation methods [START_REF] Giesl | Construction of Global Lyapunov Functions Using Radial Basis Functions[END_REF], [START_REF] Johansen | Computation of Lyapunov functions for smooth nonlinear systems using convex optimization[END_REF], graph theoretic methods [START_REF] Ban | A computational approach to Conley's decomposition theorem[END_REF], [START_REF] Kalies | An algorithmic approach to chain recurrence[END_REF], and semidefinite optimization for sum-of-squares polynomials (known as the SOS method) [START_REF] Papachristodoulou | The construction of Lyapunov functions using the sum of squares decomposition[END_REF], [START_REF] Peet | A converse sum-of-squares Lyapunov result: An existence proof based on the Picard iteration[END_REF]. † S. Hafstein is with the School of Science and Engineering at Reykjavik University, Iceland sigurdurh@ru.is. ‡ C.M. Kellett is with the School of Electrical Engineering and Computer Science, University of Newcastle, Callaghan, New South Wales, Australia Chris.Kellett@newcastle.edu.au. In this paper we consider systems described by ordinary differential equations

ẋ = f (x),x 2 R n , (1)
where we assume f (•) is twice continuously differentiable, and denote solutions to (1) by φ :

R ≥0 ⇥ R n ! R n .
One proposed approach to constructing a Lyapunov function for (1) involves partitioning the state-space into simplices (called a triangulation), defining values for the vertices of every simplex, and for every simplex taking the convex interpolation of those values. This yields a continuous and piecewise affine (CPA) function. If the values at the vertices satisfy a system-dependent set of linear inequalities, then the resulting CPA function is a Lyapunov function.

An approach using linear programming to compute feasible values for the CPA function at the simplex vertices was proposed in [START_REF] Marinosson | Lyapunov function construction for ordinary differential equations with linear programming[END_REF] with refinements in [START_REF] Hafstein | An Algorithm for Constructing Lyapunov Functions[END_REF], [START_REF] Baier | Linear programming based Lyapunov function computation for differential inclusions[END_REF], [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF]. However, since linear programming is not very efficient, the question remains if the values at the vertices cannot be fixed by more efficient means with a subsequent speedy test of the validity of the linear inequalities. In this paper, we propose using a construction from a particular converse Lyapunov theorem. Classical converse Lyapunov theorems such as those developed by Massera [START_REF] Massera | On Liapounoff's conditions of stability[END_REF] and Kurzweil [START_REF] Kurzweil | On the inversion of Ljapunov's second theorem on stability of motion[END_REF] rely on integrating solutions from the initial time to infinity. However, Yoshizawa [START_REF] Yoshizawa | Stability Theory by Liapunov's Second Method[END_REF] provided an alternate construction that involves taking the supremum over time of the norm of the solution. Initially this appears to provide no improvement towards a constructive approach, but it can be shown that this supremum is actually a maximum over a finite-time horizon. Furthermore, in many cases, this horizon may not be overly long.

We therefore propose a method for constructing CPA Lyapunov functions based on using Yoshizawa's construction for the values at the simplex vertices and subsequently verifying that this yields a true Lyapunov function by checking the validity of the linear inequalities from Theorem 2 below. In Theorem 3 below we demonstrate that this construction will always succeed if the CPA function has enough structure, i.e. the triangulation has a sufficient number of vertices, and the numerical integration of the solution trajectories delivers reliable values.

The paper is organized as follows: in Section II we describe the construction of CPA functions on a given triangulation and the linear program used to verify if a given CPA function is, in fact, a Lyapunov function. In Section III we describe the Lyapunov function construction due to Yoshizawa and describe the form of the stability estimates required. In Section IV we present three representative numerical examples and in Section V we conclude and provide some indications of future work.

II. CONTINUOUS AND PIECEWISE AFFINE LYAPUNOV FUNCTIONS

In the sequel, we will define continuous and piecewise affine (CPA) functions on suitable triangulations. For a set Ω ⇢ R n , we denote the interior of Ω by Ω • and the closure of Ω by Ω.

Definition 1: We call a finite collection

T = {S 1 , S 2 ,...,S N } of n-simplices in R n a suitable triangulation if i) S ν , S µ 2T, ν 6 = µ, intersect in a common face or not at all. ii) With D T . = [ ν S ν , D • T is a connected neighborhood of the origin. iii) If 0 2 S ν , then 0 is a vertex of S ν .
Remark 1: Property i), often called shape regularity in the theory of FEM, is needed so that we can parameterize every continuous function, affine on every simplex, by specifying its values at the vertices, cf. Remark 2. Property ii) ensures that D T is a natural domain for a Lyapunov function and, without Property iii), a function affine on each of the simplices could not have a local minimum at the origin. ⇤ For a given suitable triangulation, T , and with D T . = [ S2T S, we denote the set of all continuous functions f : D T ! R that are affine on every simplex S 2T by CPA[T ].

Remark 2: A function f 2 CPA[T ] is uniquely determined by its values at the vertices of the simplices of T .T o see this, let S ν =c o {x 0 ,x 1 ,...,x n }2T. Every point x 2 S ν can be written uniquely as a convex combination of its vertices, x = P n i=0 λ x i x i , λ x i ≥ 0 for all i =0, 1,...,n, and

P n i=0 λ x i =1 . The value of f at x is given by f (x)= P n i=0 λ x i f (x i ). Additionally, f has a representation on S ν as f (x)=w T (x -x 0 )+a for w 2 R n and a 2 R. Let rf ν . = w.
Then, as shown in [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF]Remark 9], rf ν is linear in the values of f at the vertices x 0 ,x 1 ,...,x n . ⇤ Our subsequent results will be valid on a domain D⇢R n minus a fixed arbitrarily small neighborhood of the origin. We define a CPA[T ] Lyapunov function that accounts for this.

Definition 2: Let T be a suitable triangulation and let V 2 CPA[T ] be a positive definite function. Define constants

α ⇤ 1 . =i n f x2D T \{0} V (x) |x| ,α ⇤ 2 . =s u p x2D T \{0} V (x) |x| , and r . =min x2∂D T |x|. Let ε 2 R >0 satisfy 0 <ε<rα ⇤ 1 /α ⇤ 2 (2)
and define

B ε . = {x 2 R n : |x| <ε}. If there is a constant α ⇤ 3 > 0 such that lim sup h!0+ V (x + hf (x)) -V (x) h -α ⇤ 3 |x| for all x 2 (D T \ B ε ) • we call V a CPA[T ] Lyapunov function for (1) on D T \ B ε .
The implication of a CPA[T ] Lyapunov function for (1) on D T \ B ε is slightly weaker than asymptotic stability.

Theorem 1: Given a suitable triangulation, T , and

ε 2 R >0 , assume that V : D!R ≥0 is a CPA[T ] Lyapunov function for (1) on D T \ B ε . For every c 2 R ≥0 define the sublevel set L V,c . = {x 2D T : V (x)  c} and let m . = max |x|=ε V (x) and M . =m i n x2∂D T V (x). Then, for every c 2 [m, M ) we have B ε ⇢ L V,c ⇢D •
T and there exists a T c ≥ 0 such that φ(t, L V,c ) ⇢ L V,m for all t ≥ T c .

In other words, a CPA[T ] Lyapunov function implies asymptotic stability of the set L V,m . The proof is similar to [6, Theorem 6.16] and we omit the details due to space constraints.

The following theorem and corollary provide a set of linear inequalities such that, if a given CPA function satisfies the inequalities then it is a CPA Lyapunov function.

Theorem 2: Let T be a suitable triangulation and let

V 2 CPA[T ]. Define S ν . =c o {x ν 0 ,x ν 1 ,...,x ν n }2T and let B ν 2 R >0 satisfy B ν ≥ max i,j,k=1,2,...,n x∈Sν ∂ 2 f k ∂x i ∂x j (x) .
For each S ν , for i =0, 1,...,n define the constants

E i,ν . = nB ν 2 |x i -x 0 | (|x i -x 0 | + diam(S ν )) .
Then, for every S ν such that the inequalities

0 > rV ν • f (x ν i )+|rV ν | 1 E i,ν (3)
hold for all i =0, 1,...,n, we have

0 > rV ν • f (x)
for all x 2 S ν . Corollary 1: Assume that V 2 CPA[T ] from Theorem 2 is positive definite and that the constant ε 2 R >0 fulfills [START_REF] Ban | A computational approach to Conley's decomposition theorem[END_REF]. If the inequalities [START_REF] Barbashin | Existence of smooth solutions of some linear equations with partial derivatives[END_REF] are fulfilled for all S ν 2T with

S ν \ B C ε 6 = ;, then V is a CPA Lyapunov function for (1) on D T \ B ε .
The proofs of Theorem 2 and Corollary 1 can be proved in a manner similar to [5, Theorem 1]. We omit the details.

Remark 3: The usefulness of Theorem 2 is that it reduces the verification that a function V 2 CPA[T ] is a Lyapunov function for [START_REF] Baier | Linear programming based Lyapunov function computation for differential inclusions[END_REF] to the verification of a finite number of inequalities [START_REF] Barbashin | Existence of smooth solutions of some linear equations with partial derivatives[END_REF]. Finding a candidate CPA Lyapunov function can be done as in [START_REF] Baier | Linear programming based Lyapunov function computation for differential inclusions[END_REF], [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF], [START_REF] Hafstein | An Algorithm for Constructing Lyapunov Functions[END_REF], [START_REF] Marinosson | Lyapunov function construction for ordinary differential equations with linear programming[END_REF], via linear programming. Alternatively, as in this paper, one can define V 2 CPA[T ] by computing suitable values at the vertices of the simplices of T and then verify the inequalities (3). ⇤ In the next theorem we consider CPA approximations to functions.

Definition 3: Let D⇢R n be a domain, f : D!R be a function, and T be a triangulation such that D T ⇢D . The CPA[T ] approximation g to f on D T is the function g 2 CPA[T ] defined by g(x)=f (x) for all vertices x of all simplices in T .

We additionally need that the simplices in the triangulation T are not too close to being degenerate. This property can be quantified as follows: For an n-simplex S ν . = co{x 0 ,x 1 ,...,x n }2T define its shape-matrix as X ν by writing the vectors x 1x 0 , x 2x 0 ,...,x nx 0 in its rows subsequently. The degeneracy of the simplex X ν can now be quantified through the value diam(S ν )|X -1 ν |, where |X -1 ν | is the spectral norm of the inverse of X ν (see part (ii) in the proof of [START_REF] Baier | Linear programming based Lyapunov function computation for differential inclusions[END_REF]Theorem 4.6]). Then for every R>0 there exists a δ R > 0 such that, for any triangulation T satisfying 1)

C⇢D T ⇢D, 2) max Sν 2T diam(S ν )  δ R , and 3) max Sν 2T diam(S ν )|X -1 ν |R the CPA[T ] approximation V to W on D T is a CPA Lyapunov function for (1) on D T \ B ε .
The proof of Theorem 3 follows [5, Theorem 5] and we omit the details due to space constraints.

Remark 4: Theorem 3 is more constructive than it might seem at first glance since a given triangulation T can be manipulated to deliver a new triangulation T ⇤ with smaller simplices without increasing their degeneracy. As a consequence, it is always possible to find a triangulation that admits a CPA Lyapunov function approximating a twice continuously differentiable Lyapunov function. ⇤ III. YOSHIZAWA CONSTRUCTION OF LYAPUNOV FUNCTIONS We now turn to the question of how to define the vertex values of each simplex in order to obtain a CPA Lyapunov function. We propose using a numerical approximation of a construction initially proposed by Yoshizawa in proving a converse Lyapunov theorem [START_REF] Yoshizawa | Stability Theory by Liapunov's Second Method[END_REF]. We make use of the standard function classes K 1 and KL (see [START_REF] Hahn | Stability of Motion[END_REF], [START_REF] Kellett | A compendium of comparsion function results[END_REF]).

Let the open set D⇢R n be such that the origin is contained in D. Suppose ( 1) is KL-stable on D; i.e., there exists β 2KLso that

|φ(t, x)|β(|x|,t), 8x 2D,t2 R ≥0 . (4)
It was shown in [20, Proposition 1] that KL-stability is equivalent to (local) asymptotic stability of the origin for (1) where D is contained in the basin of attraction. When D = R n , KL-stability is equivalent to global asymptotic stability of the origin for (1). We will refer to the function β 2KLof (4) as a stability estimate.

In what follows we will make use of Sontag's lemma on KL-estimates [19, Proposition 7] ([10, Lemma 7]):

Lemma 1: Given β 2K Land λ 2 R >0 , there exist α 1 ,α 2 2K 1 so that, for all s, t 2 R ≥0 α 1 (β(s, t))  α 2 (s)e -λt .
Definition 4: Given a stability estimate β 2K L , let α 1 ,α 2 2K 1 come from Lemma 1 with λ =2 . We call the function V : R n ! R ≥0 defined by

V (x) . =sup t≥0 α 1 (|φ(t, x)|)e t (5)
a Yoshizawa-Lyapunov function.

As shown in [START_REF] Teel | A smooth Lyapunov function from a class-KL estimate involving two positive semidefinite functions[END_REF] the Yoshizawa-Lyapunov function satisfies the bounds

α 1 (|x|)  V (x)  α 2 (|x|) (6)
and the decrease condition

V (φ(t, x))  V (x)e -t . (7)
The lower bound in ( 6) is trivial by considering t =0while the upper bound follows from

V (x)  sup t≥0 α 1 (β(|x|,t))e t  α 2 (|x|)e -2t+t  α 2 (|x|). (8)
As shown in [START_REF] Teel | A smooth Lyapunov function from a class-KL estimate involving two positive semidefinite functions[END_REF], when the vector field of ( 1) is locally Lipschitz the function V (•) is continuous (except possibly at the origin) and there exists a time T (x) so that

V (x) = max t2[0,T (x)] α 1 (|φ(t, x)|)e t . (9)
We can calculate the time T (x) explicitly as [20, Claim 2]

T (x)=-ln ✓ V (x) α 2 (|x|) ◆ +1,x 6 =0 (10)
and with the upper and lower bounds on V we see that

0  T (x) -ln ✓ α 1 (|x|) α 2 (|x|) ◆ +1 =l n ✓ α 2 (|x|) α 1 (|x|) ◆ +1. (11)
Recall that our intention is to calculate V (x) for each x that is a simplex vertex. In order to do this, we clearly need a solution to (1) from each such x. As a closed form solution is generally not available, we will resort to numerical integration in order to calculate V (x) given by [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF]. For this approach to be numerically tractable, we require that the time horizon T (x) given by ( 11) not be too large. We present two examples of stability estimates and derive T (x) in each case.

A. Example 1 -Exponentially Stable Estimates

Suppose the stability estimate is given by |φ(t, x)|α(|x|)e -µt ,µ > 0, where α(s) ≥ s. Then, for all s 2 R ≥0 , we can define α 1 (s) . = s 2/µ , and α 2 (s) . =(α(s)) 2/µ so that, for all x 2 R n and t 2 R ≥0 ,

α 1 (α(|x|)e -µt )  (α(|x|)) 2/µ e -2t = α 2 (|x|)e -2t .
We therefore see that, in this case, an upper bound for the time horizon to optimize over is given by

T (x)  2 µ ln ✓ α(|x|) |x| ◆ +1 (12)
where the assumption that α(s) ≥ s for all s 2 R ≥0 guarantees that T (x) ≥ 1. If α(s)=Ms for some M>1, then an upper bound for the time horizon necessary to optimize over is independent of the point x and is given by

T (x)=T = -lnM -2/µ +1= 2
µ lnM +1. In Section IV-A this bound will be used to compute values of the Yoshizawa-Lyapunov function for a linear system.

B. Example 2

With the functions α 1 ,α 2 2K 1 given by

α -1 1 (s) . = e s -1,α 2 (s)=Ms we capture functions β 2KLsatisfying β(s, t)  exp(Mse -2t ) -1 (13)
and the optimization horizon bound is given by

T (x)  ln ✓ M |x| ln(1 + |x|) ◆ +1.
The horizon length grows with increasing |x| but not too quickly. For example, with M = 10: |x| =1yields T (x)= 3.67 and |x| = 100 yields T (x)=6.38. Remark 5: There are two difficulties we encounter in trying to calculate [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF]. The first difficulty lies with finding a stability estimate β 2K Lor even with verifying that a particular stability estimate such as [START_REF] Lyapunov | The general problem of the stability of motion[END_REF] holds for a particular system (1). There seems to be little that can be done to circumvent this problem.

The second difficulty is that Sontag's lemma on KLestimates is not constructive and, to the best of the authors' knowledge, given an arbitrary β 2K L , there are currently no constructive techniques for finding α 1 ,α 2 2K 1 .

Remark 6: As previously mentioned, in [START_REF] Teel | A smooth Lyapunov function from a class-KL estimate involving two positive semidefinite functions[END_REF] it was demonstrated that the Yoshizawa-Lyapunov function ( 5) is continuous except possibly at the origin. However, in order to make use of Theorem 3 we would clearly prefer that the Yoshizawa-Lyapunov function be twice continuously differentiable. Standard smoothing techniques can be applied to the Yoshizawa-Lyapunov function to obtain a smooth function, however these techniques do not generally provide us with an easily calculable function such as [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF]. Rather, we would prefer to demonstrate directly that the Yoshizawa-Lyapunov function inherits the regularity of the vector field of (1). It is known that, in general, there exists a Lyapunov function that inherits the regularity property of the vector field defining (1) (see [START_REF] Barbashin | Existence of smooth solutions of some linear equations with partial derivatives[END_REF]), however, it remains to be shown that that the Yoshizawa-Lyapunov function (5) inherits the regularity property of the vector field. This is a subject of future work.

IV. NUMERICAL EXAMPLES

While the Yoshizawa construction (9) provides an exact value for a Lyapunov function at each point, this assumes we have an exact solution of the differential equation. In this section we present three numerical examples of CPA Lyapunov functions. In each case, we first define a triangulation on a region of the state space that includes the origin in its interior. For technical reasons, we excise a small neighborhood of the origin.

For each example we calculate a stability estimate and, with the triangulation defined, we then use a simple Euler integration scheme to approximate the values of the Yoshizawa-Lyapunov function ( 9) at all vertices of all simplices. We then check the inequalities (3) to verify that the function defined by taking the convex interpolation on each simplex of the Yoshizawa-Lyapunov function values is, in fact, a CPA-Lyapunov function.

A. Example 3 -Linear System

Consider the linear system

ẋ = Ax =  11 -5 -3 . ( 14
)
We observe that the origin is globally exponentially stable as the eigenvalues are at -1 ± i and, by solving the Lyapunov equation A T P + PA = -Id, a Lyapunov function is given by

V (x)=x T Px = x T  4.51 10 .5 x. (15)
By explicitly calculating the solutions of [START_REF] Marinosson | Lyapunov function construction for ordinary differential equations with linear programming[END_REF] we see that the system satisfies the stability estimate

|φ(t, x)|7|x|e -t , 8x 2 R 2 ,t2 R ≥0 . (16)
From Section III-A, with α(s)=7s and µ =1we see that α 1 (s)=s 2 , α 2 (s) = 49s 2 , and T (x)=T =5. According to the above proposed procedure, we define a triangulation and define the values at the simplex vertices by approximating [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF] to obtain a continuous and piecewise affine function V 1 (x) for ( 14) as shown in Figure 1. The neighborhood of the origin {x 2 R 2 : |x| < 0.048}, is excluded. It is straightforward to numerically verify the inequalities (3) to conclude that the obtained function, V 1 (x), is in fact a CPA-Lyapunov function.

The function V (x) given by ( 15) has a similar though slightly different shape. Level curves for V 1 (x) are shown in Figure 2 and level curves for V (x) are shown in Figure 3 for comparison.

B. Example 4 -Simple Nonlinear System

Consider the system

ẋ = -x 3 (17)
which has solution We observe that the norm of the solution is in fact a KL function and, consequently, immediately provides a stability estimate. We can verify that the functions

φ(t, x)= x p 1+2x 2 t , 8x 2 R,t 2 R ≥0 . (18)
α 1 (s)=α 2 (s) . = ⇢ 0 ,s =0 s exp -1 s 2 ,s > 0 (19) are such that α 1 ✓ s p 1+2s 2 t ◆  α 2 (s)e -2t . (20)
As before, we define a triangulation and calculate the values at the vertices by approximating [START_REF] Giesl | Revised CPA method to compute Lyapunov functions for nonlinear systems[END_REF]. A convex interpolation of these values on each simplex then yields a CPA function and the inequalities (3) are used to verify that the calculated CPA function is a CPA-Lyapunov function.

We note that, for any p 2 Z ≥1 and c 2 R >0 , a Lyapunov function for ( 17) is given by

V (x)=cx 2p , 8x 2 R. ( 21
)
Figure 4 shows the CPA-Lyapunov function V 2 (x) for system ( 17

C. Example 5 -Nonlinear System

Consider the two-dimensional nonlinear system given by

ẋ1 = -x 2 -(1 -x 2 1 -x 2 )x 1 ẋ2 = x 1 -(1 -x 2 1 -x 2
2 )x 2 .

(22

)
This system has the unit circle as a periodic orbit and the origin as a locally asymptotically stable equilibrium. On any compact subset of the unit ball, the simple quadratic V (x) . and T (x)=T =1.

For this example,with R =0 .94478 and using the numerical procedure previously outlined, a CPA-Lyapunov function V 3 (x) of system (22) was computed and is shown in

V. C ONCLUSIONS

In this paper we have presented a novel technique for the numerical construction of Lyapunov functions given a stability estimate in the form of a KL-bound on the norm of system trajectories. For a suitable triangulation of the statespace, at each simplex vertex we calculate the value of a Lyapunov function construction due to Yoshizawa [START_REF] Yoshizawa | Stability Theory by Liapunov's Second Method[END_REF]. From these values, we then define a CPA function on the domain minus an arbitrarily small neighborhood of the origin. We can verify that the CPA function thus defined is a Lyapunov function (Corollary 1). If the CPA function thus defined is not a Lyapunov function (i.e., does not satisfy (3)), we can refine the triangulation and construct a new CPA function. If the Yoshizawa-Lyapunov function is twice continuously differentiable, then this process of refining the triangulation then yields a CPA Lyapunov function in a finite number of steps (Theorem 3).

Theorem 3 :

 3 Let C, D⇢R n be simply connected compact neighborhoods of the origin such that C • = C, D • = D, and C⇢D • . Assume that W 2 C 2 (D) is a Lyapunov function for (1). Set r . =min x2∂C |x| and let ε 2 R >0 satisfy ε<r• inf x2D\{0} W (x)/|x| sup x2D\{0} W (x)/|x| .

Fig. 1 .

 1 Fig. 1. CPA Lyapunov function V 1 (x) for system (14).

Fig. 2 .Fig. 3 .

 23 Fig. 2. Level curves of V 1 (x) for values 0.189, 0.378, 0.567, 0.756, and 0.945.

  Figure 4 shows the CPA-Lyapunov function V 2 (x) for system (17) for 4.8 ≥| x|≥0.04332. For comparison, Figure 4 also shows the known Lyapunov function (21) with p =2, c =0.01.

Fig. 4 .

 4 Fig.4. Lyapunov functions V (x)=0.01x 4 (green curve) and V 2 (x) (red curve) for system[START_REF] Peet | A converse sum-of-squares Lyapunov result: An existence proof based on the Picard iteration[END_REF].

  is a known Lyapunov function.Fix R 2 (0, 1). Then, for any initial conditions satisfyingx 2 1 + x 2 2  Rwe have the stability estimate|φ(t, x)||x|e -(1-R)t(24)and, from Section III-A, we can calculate α 1 (s)=α 2 (s)=s 2/(1-R)

  Figure 5 for |x|2[0.012,R]. For comparison, Figure 5 also shows the known Lyapunov function (23).

Fig. 5 .

 5 Fig. 5. Difference between Lyapunov functions V (x) (upper arrow points) and V 3 (x) (lower arrow points) for system (22).

  Kellett is supported by ARC Future Fellowship FT1101000746 and by the Alexander von Humboldt Foundation.

* H. Li is with the Mathematical Institute, University of Bayreuth, 95440 Bayreuth, Germany Huijuan.Li@uni-bayreuth.de. Li is supported by the EU Initial Training Network "Sensitivity Analysis for Deterministic Controller Design -SADCO".


Thermoelasticity of non-simple oriented materials

Czeslaw Wozniak

To cite this version: TaKxce OT rpa,qnefrToB TehfnepaTypbr. noHKTsie opaefuHposaHHor0 MaTep~anao6oioti~eHoHaTep~~Kuenpo6neM6InyTeMBBeAeHHR,He3aBHC~OOTAO6aBOSHbIXCTCneHC~ cao6onar,Ta1cxe~o6aBow%r~ nonet,oTo6pamaHx.twiXpacnpeAeneH~e TeMrrepaTyp.

  Czeslaw Wozniak. Thermoelasticity of non-simple oriented materials. International Journal of Engineering Science, 1967, 5 (8), pp.605-612. hal-00944400


Gustav Nordh

Bruno Zanuttini

Frozen Boolean partial co-clones

We introduce and investigate the concept of frozen partial co-clones. Our main motivation for studying frozen partial co-clones is that they have important applications in complexity analysis of constraints. The frozen partial coclones lie between the co-clones and partial co-clones in the sense that the partial co-clone lattice is a refinement of the frozen partial co-clone lattice, which in turn is a refinement of the co-clone lattice. We concentrate on the Boolean domain and determine large parts of the frozen partial coclone lattice.

Introduction

A clone is a composition closed set of operations containing all projections. An operation f preserves a relation R (or is a polymorphism of R) if f applied coordinate-wise to any tuples from R gives a tuple in R. Given a set of relations Γ, P ol(Γ) is the set of operations preserving all relations in Γ. P ol(Γ) is always a clone and any clone can be defined as P ol(Γ) for a set of relations Γ. Inv(F ) is the set of relations that are preserved by (invariant under) all the operations in F .

An n-ary relation R has a primitive positive (p.p.) definition (also called an implementation) in a set of relations Γ if R is the set of models of an existentially quantified conjunction of atomic formulas over Γ ∪ {=} (also called a Γ-formula), i.e., R(x 1 , . . . , x n ) ≡ ∃X i R i (x i1 , . . . , x it ) where each R i ∈ Γ ∪ {=}. Sets of relations closed under p.p. definability are called co-clones and the least co-clone containing Γ is denoted by Γ . There is a Galois connection between clones and sets of relations closed under p.p. definability (i.e., co-clones). In particular, given two sets of relations Γ 1 and Γ 2 , then Γ 1 ⊆ Γ 2 if and only if P ol(Γ 2 ) ⊆ P ol(Γ 1 ), and for any set of relations Γ we have Γ = Inv(P ol(Γ)). For more information on clones, co-clones, and the aforementioned Galois connection, we refer the reader to the books [START_REF] Lau | Function Algebras on Finite Sets[END_REF][START_REF] Pöschel | Funktionen-und Relationenalgebren[END_REF].

The CSP(Γ) problem, where Γ is a set of relations (also called a constraint language), is the problem of deciding if a given set of variables subject to a set of constraints (given by atomic formulas over Γ) is satisfiable. The problem of classifying the computational complexity of CSP(Γ) with respect to Γ is an important open problem. The (so far) most successful approach to attack this problem is an algebraic approach which heavily relies on the following result.

Theorem 1 ( [START_REF] Jeavons | On the algebraic structure of combinatorial problems[END_REF][START_REF] Jeavons | Closure properties of constraints[END_REF])

If Γ 1 is finite and Γ 1 ⊆ Γ 2 , then CSP(Γ 1 ) is polynomial-time reducible to CSP(Γ 2 ).
Note that because of the Galois connection between clones and co-clones this result can be reformulated into: If Γ 1 is finite and P ol(Γ 2 ) ⊆ P ol(Γ 1 ), then CSP(Γ 1 ) is polynomial-time reducible to CSP(Γ 2 ). To illustrate the power of this result, note that Schaefer's [START_REF] Schaefer | The complexity of satisfiability problems[END_REF] seminal complexity classification (separating the cases in P from the NP-complete cases) for CSP(Γ) over the Boolean domain follows trivially from this result and Post's classification of Boolean clones [START_REF] Post | The two-valued iterative systems of mathematical logic[END_REF] (i.e., Post's lattice). For more information on the connection between clones, co-clones, and the CSP(Γ) problem, we refer the reader to the survey articles [START_REF] Böhler | Playing with Boolean blocks, part I: Post's lattice with applications to complexity theory[END_REF][START_REF] Böhler | Playing with Boolean blocks, part II: Constraint satisfaction problems[END_REF][START_REF] Krokhin | Functions of multiple-valued logic and the complexity of constraint satisfaction: A short survey[END_REF].

This result is not so useful for more fine grained complexity analysis of CSP(Γ) (and similar) problems. The reason is that it does not preserve the size of the problem instances (in terms of the number of variables). The crux is that in order for the proof of Theorem 1 to work, the existential quantifiers in p.p. definitions are eliminated by introducing new variables. So the reduction in Theorem 1 maps a conjunction of m constraints on n variables to a conjunction of bm constraints on n + am variables, where a and b are constants which depend on the languages Γ 1 and Γ 2 .

As an example of a problem for which this is too coarse a reduction, compare 3-SAT (i.e., CSP(Γ 3SAT ) where Γ 3SAT consists of the relations corresponding to clauses on at most three variables) with 1-in-3-SAT (i.e., CSP(Γ 1/3 ) where Γ 1/3 consists of the relation {001, 010, 100}). By Post's lattice it is easy to verify that Γ 1/3 = Γ 3SAT and hence, CSP(Γ 3SAT ) is polynomial-time equivalent to CSP(Γ 1/3 ) according to Theorem 1. Despite this, 3-SAT (solvable in time O(1.473 n ) [START_REF] Brueggemann | An improved deterministic local search algorithm for 3-SAT[END_REF][START_REF] Dantsin | A deterministic (2-2/(k+1)) n algorithm for k-SAT based on local search[END_REF]) seems to be a much harder problem than 1-in-3-SAT which can be solved in time O(1.1003 n ) [START_REF] Byskov | New algorithms for exact satisfiability[END_REF] (where n is the number of variables). Hence, it is clear that to get a better understanding of the complexity of CSP(Γ) (and similar) problems, we need reductions/implementations/tools where the blow-up in instance size can be controlled.

With this in mind we consider partial co-clones instead. A partial clone is a composition closed set of partial functions containing all (total) projections. A partial operation f preserves a relation R (or is a partial polymorphism of R) if f applied coordinate-wise to any tuples from R gives a tuple in R whenever f is defined on all the arguments. pP ol(Γ) is the set of (partial) operations preserving all relations in Γ. Partial co-clones can be defined as the sets of relations that are closed under p.p. definitions not using existential quantification, and the least partial co-clone containing Γ is denoted by Γ p . There is also a Galois connection between partial co-clones and partial polymorphisms, in particular we have the following result.

Theorem 2 ([9, 8, 17]) Let Γ 1 and Γ 2 be sets of relations. Then Γ 1 p ⊆ Γ 2 p if and only if pP ol(Γ 2 ) ⊆ pP ol(Γ 1 ).

For formal definitions and more information on partial clones and partial co-clones we refer the reader to [START_REF] Lau | Function Algebras on Finite Sets[END_REF].

Since partial co-clones do not utilize existential quantification, it can be verified that if we replace co-clones by partial co-clones in Theorem 1 then the reductions are sizepreserving (in terms of the number of variables). Hence, if

Γ 1 p ⊆ Γ 2 p and Γ 1 is finite then CSP(Γ 1 ) is solvable in time O(f (n)) if CSP(Γ 2 ) is solvable in time O(f (n))
(where n is the number of variables). This result seems more useful than what it is since the structure of partial coclones is very complicated and only small portions of the lattice are known even for the Boolean domain [START_REF] Haddad | On intervals of partial clones of Boolean partial functions[END_REF]. We remark in passing that there are other applications of partial clones in complexity analysis of CSP(Γ) problems [START_REF] Schnoor | New algebraic tools for constraint satisfaction[END_REF].

What we would like to have is a concept that combines the good features of partial co-clones (e.g., size-preserving reductions for CSP(Γ)) and the good features of co-clones (e.g., simpler structure especially over the Boolean domain). In this paper we propose and investigate such a concept that we call frozen partial co-clones. Frozen partial co-clones can be defined as the sets of relations that are closed under p.p. definitions where only the variables that take the same value in every model of the p.p. definition/formula (so called frozen variables) can be existentially quantified. The least frozen partial co-clone containing Γ is denoted Γ f r . Hence, the frozen partial co-clones lie between the co-clones and partial co-clones in the sense that the partial co-clone lattice is a refinement of the frozen partial co-clone lattice, which in turn is a refinement of the co-clone lattice. Moreover, if we replace co-clones by frozen partial co-clones in Theorem 1 then it can be verified that if

Γ 1 f r ⊆ Γ 2 f r and Γ 1 is finite then CSP(Γ 1 ) is solv- able in time O(f (n + |D|)) if CSP(Γ 2 ) is solvable in time O(f (n)) (
where D is the domain). The point is that all the variables that are frozen to the same domain element can be replaced by a single variable, and hence, at most |D| extra variables need to be introduced when eliminating the existential quantifiers.

In this paper we focus exclusively on the Boolean domain. As expected, there is a connection between frozen partial co-clones and partial polymorphisms that we sketch (for the Boolean domain) in Section 2. In Sections 3 and 4 we determine large portions of the frozen (Boolean) partial co-clone lattice which indeed is significantly simpler than the partial co-clone lattice.

Frozen partial co-clones and partial polymorphisms

In this section we give some definitions and preliminary remarks before sketching the connection between frozen (Boolean) partial co-clones and frozen (Boolean) partial polymorphisms. If ϕ is a formula, then V ars(ϕ) denotes the set of variables occurring in it, and M(ϕ) denotes the set of all assignments to V ars(ϕ) which satisfy ϕ (i.e., the models of ϕ). The relations {0} and {1} are denoted by F and T , respectively. Atomic formulas are sometimes written in prefix notation (e.g., R(x 1 , . . . , x n )), or infix notation (e.g., x 1 = x 2 ), depending on the context. Given a function f , dom(f ) denotes the domain of f (i.e., the set of tuples t i for which f (t i ) is defined), and f is a subfunction of g if dom(f ) ⊆ dom(g) and f (t i ) = g(t i ) for all t i ∈ dom(f ).

Definition 3 (frozen variable) Let ϕ be a formula and let x ∈ V ars(ϕ). Then x is said to be frozen in ϕ if ϕ |= T (x) or ϕ |= F (x). In other words, x is frozen in ϕ if it is assigned the same value by all its models. Definition 4 (frozen implementation) Let Γ be a set of relations and R an n-ary relation. Then Γ freezingly implement R if there is a p.p. definition R(x 1 , . . . , x n ) ≡ ∃Xϕ such that ϕ is a conjunction of atomic formulas over Γ ∪ {=}, V ars(ϕ) ⊆ X ∪{x 1 , . . . , x n }, and every variable in X is frozen in ϕ.

Note that frozen implementations are slightly less general than so-called faithful implementations [5, page 34]. But the latter are not suitable for our purposes since they blow up instance sizes.

Definition 5 (frozen partial co-clone) Let Γ be a set of relations. The frozen partial co-clone generated by Γ, written It is easy to see that frozen implementations compose together, in the sense that if Γ freezingly implements every relation in Γ ′ and Γ ′ freezingly implements R, then Γ freezingly implements R. It is also easy to see that frozen partial co-clones ordered by set inclusion form a lattice.

Example 6

The relation R ′ = {00, 10} is in {R, T } f r with R = {000, 001, 110} and T = {1}. This is because

R ′ (x 1 , x 2 ) ≡ ∃x 3 R(x 2 , x 2 , x 3 ) ∧ R(x 2 , x 2 , x 1 ) ∧ T (x 3 )
where x 3 is frozen to 1.

As another example, consider the relations R 1 = {01, 10} and R 2 = {0100, 0110, 1000, 1111}. Then

R ′ = {010, 011, 100} is in {R 1 , R 2 } f r , as shown by R ′ (x 1 , x 2 , x 3 ) ≡ ∃x 4 R 1 (x 1 , x 2 ) ∧ R 2 (x 1 , x 2 , x 3 , x 4 ) (x 4
is frozen to 0 by the conjunction).

Definition 7 (determined) If Γ is a set of relations such that there is a Γ-formula ϕ in which x is frozen to d ∈ {0, 1}, then we say that d is determined in Γ.

Proposition 8 Let Γ be a set of relations. Then d ∈ {0, 1} is determined in Γ if and only if {d} ∈ Γ f r .

Proof: Without loss of generality assume d = 1. If 1 is determined in Γ, then there is a Γ-formula ϕ and a variable x T ∈ V ars(ϕ) such that ∃Xϕ |= T (x T ) and thus ϕ |= T (x T ). Let m be a model of ϕ with a maximum number of variables being assigned 1. Identify all variables in ϕ that are assigned 1 by m to x T , resulting in ϕ ′ . If V ars(ϕ ′ ) = {x T }, then T (x T ) ≡ ϕ ′ , and thus T ∈ Γ f r . Otherwise, there is a variable x F ∈ V ars(ϕ ′ ) \ {x T }. Identify all variables in V ars(ϕ ′ ) \ {x T } to x F , resulting in ϕ ′′ . Then, T (x T ) ≡ ∃x F ϕ ′′ , and T ∈ Γ f r since x F is frozen in ϕ ′′ .

The converse follows directly from the definition of determined constants. ✷ Definition 9 (frozen partial polymorphisms) A k-ary (partial) function f ∈ pP ol(Γ) is said to be frozen if it is defined on every all-d-tuple (of length k) for which d is determined in Γ and f (d, d, . . . , d) = d. We define f rP ol(Γ) = {f ∈ pP ol(Γ) | f is f rozen}, and say that f rP ol(Γ) are the frozen polymorphisms of Γ.

Note that if no d is determined in Γ, then pP ol(Γ) = f rP ol(Γ) and Γ p = Γ f r .

Lemma 10 Let Γ be a frozen partial co-clone, then any kary (partial) function f ∈ pP ol(Γ) is a subfunction of a k-ary (partial) function g ∈ f rP ol(Γ).

Proof: It is sufficient to prove that if Γ is a frozen partial co-clone such that d ∈ {0, 1} is determined in Γ, then any f ∈ pP ol(Γ) is a subfunction of a g ∈ pP ol(Γ) such that d = (d, d, . . . , d) ∈ dom(g) and g(d, d, . . . , d) = d. Note that there can be no f ∈ pP ol(Γ) such that f (d, . . . , d) = d since this contradicts the fact (observed in Proposition 8) that {d} ∈ Γ.

Let f ∈ pP ol(Γ) be a k-ary function with dom(f ) = {t 1 , . . . , t j } such that the k-tuple d = (d, . . . , d) is not in dom(f ). We define the function f d with dom(f d ) = dom(f ) ∪ {d} such that f d (d) = d with the goal of showing that f d ∈ pP ol(Γ). Assume to the contrary that there is a relation R d ∈ Γ such that R d is not preserved by f d . This means that there are k

(not necessarily distinct) tuples t 1 , . . . , t k ∈ R d such that (f d (t 1 [1], . . . , t k [1]), . . . , f d (t 1 [m], . . . , t k [m])) /
∈ R d where m is the arity of R d . By the definition of f d and the fact that f ∈ pP ol(Γ) we know that at least one of the tuples (t 1 [j], . . . , t k [j]) = d. We can without loss of generality assume an R d such that (t

1 [1], . . . , t k [1]) = d and all other (t 1 [j], . . . , t k [j]) = d.
Since {d} ∈ Γ we have that R(y 1 , . . . , y m-1 ) ≡ ∃xR d (x, y 1 , . . . , y m-1 ) ∧ (x = d) is in Γ and thus is preserved by f . As a consequence, we have Proof:

(f (t 1 [2], . . . , t k [2]), . . . , f (t 1 [m], . . . , t k [m])) ∈ R and hence, (d, f (t 1 [2], . . . , t k [2]), . . . , f (t 1 [m], . . . , t k [m])) = (f d (t 1 [1], . . . , t k [1]), . . . , f d (t 1 [m], . . . , t k [m])) ∈ R
Since Γ ⊆ Γ f r we obviously have f rP ol( Γ f r ) ⊆ f rP ol(Γ). For the other direction, assume towards contradiction that f ∈ f rP ol(Γ) but f / ∈ f rP ol( Γ f r ). Then there is a relation R ∈ Γ f r such that f / ∈ f rP ol(R) and a Γ-formula ϕ such that R ≡ ∃Xϕ where the variables X are frozen in ϕ. This is a contradiction since f ∈ f rP ol(Γ) implies that f ∈ f rP ol(M(ϕ)) and hence f ∈ f rP ol(R) (since f is frozen). ✷ Theorem 12 Let Γ 1 and Γ 2 be sets of relations, then

Γ 1 f r ⊆ Γ 2 f r if and only if f rP ol(Γ 2 ) ⊆ f rP ol(Γ 1 ).
Proof: If there is a (partial) function f such that f ∈ f rP ol(Γ 2 ) and f / ∈ f rP ol(Γ 1 ), then assume (with the aim of reaching a contradiction) that Γ 1 f r ⊆ Γ 2 f r . From the fact that f ∈ f rP ol(Γ 2 ) and f / ∈ f rP ol(Γ 1 ) we get that there is a relation R ∈ Γ 1 such that R / ∈ Γ 2 and f / ∈ pP ol(R). By assumption Γ 1 f r ⊆ Γ 2 f r , and hence there is a frozen implementation of R using the relations in Γ 2 . So, there is a Γ 2 -formula ϕ such that R ≡ ∃Xϕ where the variables in X are frozen in ϕ. But ϕ is a Γ 2 -formula and f ∈ f rP ol(Γ 2 ), so f ∈ pP ol(M(ϕ)). Moreover, as a frozen polymorphism of Γ 2 , f is defined on all columns of M(ϕ) corresponding to variables X (Proposition 8), and finally f ∈ pP ol(R). Now we continue with the other direction. If Γ 1 f r ⊆ Γ 2 f r , then since obviously Γ f r p = Γ f r for any language Γ, we have Γ 1 f r p ⊆ Γ 2 f r p . Hence, by Theorem 2 there exists a (partial) function f ′ such that f ′ ∈ pP ol( Γ 2 f r ) and f ′ / ∈ pP ol( Γ 1 f r ). By Lemma 10, we have that f ′ is a subfunction of a (partial) function

f ∈ f rP ol( Γ 2 f r ). It is clear that f / ∈ f rP ol( Γ 1 f r ) since not even the subfunction f ′ is in pP ol( Γ 1 f r ).
Hence, by Proposition 11, we get f ∈ f rP ol(Γ 2 ) and f / ∈ f rP ol(Γ 1 ). ✷

Corollary 13 Let Γ be a set of relations. Then Inv(f rP ol(Γ)) = Γ f r .

Co-clones covered by a single frozen partial co-clone

In this section we study Boolean co-clones C such that for any set of relations Γ with Γ = C we have C = Γ f r . We say that such a co-clone is covered by a single frozen partial co-clone. We show that a large number of co-clones are covered by a single frozen partial co-clone, and hence a large part of the lattice of frozen partial co-clones is identical to the corresponding part of the lattice of co-clones. The lattice of Boolean co-clones is visualised in Figure 1, where the co-clones colored grey are covered by a single frozen partial co-clone. For explanations of the notation1 used in the lattice, we refer the reader to [START_REF] Böhler | Playing with Boolean blocks, part I: Post's lattice with applications to complexity theory[END_REF][START_REF] Böhler | Playing with Boolean blocks, part II: Constraint satisfaction problems[END_REF].

We define the most relevant co-clones here. A majority operation is a ternary operation maj satisfying maj(x, x, y) = maj(x, y, x) = maj(y, x, x) = x for all x, y ∈ {0, 1}. Similarly a minority operation is a ternary operation minor satisfying minor(x, x, y) = minor(x, y, x) = minor(y, x, x) = y for all x, y ∈ {0, 1}. The binary operations max and min return the maximum and minimum of their arguments, respectively. The co-clones ID 2 , ID 1 , and IM 2 are defined as follows: ID 2 = Inv({maj}), ID 1 = Inv({maj, minor}), and IM 2 = Inv({max, min}).

In [START_REF] Haddad | On intervals of partial clones of Boolean partial functions[END_REF] it is proved that there are 25 partial co-clones2 pC such that pC ⊆ IM 2 = Inv({max, min}), and there are 33 partial co-clones pC such that pC ⊆ ID 1 = Inv({maj, minor}). Since all co-clones C such that C ⊆ IM 2 or C ⊆ ID 1 are covered by a single frozen partial coclone, there are only 8 frozen partial co-clones f C ⊆ IM 2 and 6 frozen partial co-clones f C ⊆ ID 1 . This suggests that the lattice of frozen partial co-clones is significantly less complex than the partial co-clone lattice. Nevertheless, as expected the frozen partial co-clone lattice is more complex than the ordinary co-clone lattice. In particular we show in Section 4 that the co-clone ID 2 splits into 13 frozen partial co-clones. Moreover, it seems that none of the white co-clones in Figure 1 is covered by a single frozen partial co-clone.

The covering proofs make heavy use of the results in [START_REF] Creignou | Structure identification of Boolean relations and plain bases for co-clones[END_REF] which for every Boolean co-clone C gives a set of relations Γ such that Γ p = C. In particular, it is shown in [START_REF] Creignou | Structure identification of Boolean relations and plain bases for co-clones[END_REF] that ID 1 = { =, T, F } p . The covering proofs are numerous so due to space constraints we can only present a single illustrative case.

  d and we have a contradiction. Thus, f d ∈ pP ol(Γ). ✷ Proposition 11 Given a set of relations Γ, then f rP ol(Γ) = f rP ol( Γ f r )

Figure 1 .

 1 Figure 1. The lattice of Boolean co-clones.

Proposition 14 ID 1 =

 1 Inv({maj, minor}) = { =, T, F } p is covered by a single frozen partial coclone. Proof: Let Γ ⊆ ID 1 with Γ ⊆ ID and Γ ⊆ IR 2 (i.e., Γ = ID 1 ), where ID = = and IR 2 = {F, T } . First note that since Γ = ID 1 = { =, T, F } , we trivially have that F and T are determined in Γ. Thus, according to Proposition 8 we have {T, F } ⊆ Γ f r . Now, since Γ ⊆ IR 2 there is a relation R ∈ Γ such that

The notation used is different from Post's notation, but is now standard in the Boolean CSP area.

The results in[START_REF] Haddad | On intervals of partial clones of Boolean partial functions[END_REF] are presented in terms of partial clones.

The mnemonics are: subscripts represent the arities of the relations in the basis, and superscripts represent the nature of these relations, in the same order (p stands for positive, etc.).

Acknowledgements

We thank Henning and Ilka Schnoor and Magnus Wahlström for interesting discussions on the topic of this paper, and Steffen Reith for providing the visualisation of the lattice of Boolean co-clones in Figure 1.

Supported by ANR grant CANAR (ANR-06-BLAN-0383-02).

R /

∈ IR 2 . Then there is ϕ ≡ R of the form:

Assume no j and no k is in I, which can be ensured by propagating unary constraints. From R / ∈ IR 2 we know that R is nonempty and I = ∅. Let m ∈ R, that is, a model of ϕ, and let P = {x i | i ∈ I, m(x i ) = 1} ∪ {y i | i ∈ I, m(y i ) = 1} and N = {x i | i ∈ I, m(x i ) = 0} ∪ {y i | i ∈ I, m(y i ) = 0}. Because I = ∅ we have at least one { =}-constraint and so, P, N = ∅. Moreover, obviously every { =}-constraint in ϕ is between a variable in P and one in N . Now identify all the variables in P to a single variable p and all those in N to n. Then clearly the resulting formula is logically equivalent to (p = n) ∧ j∈J F (x j ) ∧ k∈K T (x k ), and we get a frozen implementation of = by existentially quantifying over every x j and x k (j ∈ J, k ∈ K).

Finally, { =, T, F } ⊆ Γ f r and so, Γ f r = ID 1 . ✷ Theorem 15 Each co-clone colored grey in Figure 1 is covered by a single frozen partial co-clone.

Structure of ID 2

We begin by introducing the basic relations and the 13 frozen partial co-clones in ID 2 = Inv({maj}) (i.e., the frozen partial co-clones Γ f r such that Γ f r ⊆ ID 2 and Γ = ID 2 ). We then prove that these 13 frozen partial coclones cover ID 2 (i.e., Γ f r equals one of these 13 frozen partial co-clones for any Γ such that Γ = ID 2 ). Finally, we prove that these 13 frozen partial co-clones are all distinct. We remark that the lattice of partial co-clones in ID 2

has not yet been classified [START_REF] Haddad | On intervals of partial clones of Boolean partial functions[END_REF]. Hence, the results in this section can also be seen as a step towards such a classification.

Definition 16 (relations in ID 2 ) R p 2 is the relation defined by (x 1 ∨ x 2 ) (2 stands for binary and p for positive). Similarly, R n 2 is the relation defined by (¬x 1 ∨ ¬x 2 ), R i 2 is the relation defined by (¬x 1 ∨x 2 ) (i stands for implicative), and R = 2 is the relation defined by

Definition 17 [frozen partial co-clones in ID 2 ] We define the following frozen partial co-clones:

The inclusion structure among these frozen partial coclones is given in Figure 2. Most of the inclusions are obvious. The main difficulty is to show that they cover all of ID 2 . Due to space constraints and the fact that all the covering proofs are quite similar, we only present one of them. Note that the proof is similar in spirit to the constructions in [START_REF] Schnoor | New algebraic tools for constraint satisfaction[END_REF].

Proof: Given a relation R ≡ M(ϕ), then R |{xi,...,xj } denotes the projection of R onto the coordinates corresponding to the variables {x i , . . . , x j } in ϕ.

From R ∈ Γ p 4 it follows that R has a p.p. definition of the form ∃X, i∈I R p 4 (x i1 , . . . , x i4 ) where the variables in X are frozen. Write R i for R |{xi1,...,xi4} and R p 4 (X i ) for R p 4 (x i1 , . . . , x i4 ). We claim that there is an i 0 ∈ I such that R i0 = R p 4 (X i0 ) and all x i0j 's are pairwise different. Assume to the contrary that for all i, R i = R p 4 (X i ). By construction it follows R i ⊂ R p 4 (X i ) (since at least the {R p 4 }-constraint acts on x i1 , . . . , x i4 ). But a case study on the tuples in R p 4 (X i ) \ R i shows that this entails R i ∈ ID 1 . Since from the definition of R it follows R ≡ ∃X, i∈I R i , we get R ∈ ID 1 , a contradiction. Now consider the relation obtained from R by applying the following transformations maximally while preserving R i0 = R p 4 (X i0 ) = {0110, 1001, 1100}:

∈ {x i01 , . . . , x i04 } to 0 or 1 (using F or T with Proposition 8) and existentially quantify over it, e.g., if R |{xi 0 1,...,xi 0 4,x ′ } = {01100, 01101, 10010, 11000}, freeze x ′ to 0 and existentially quantify over it.

When none of these operations can be applied any more, it is easily verified that all remaining x ′ / ∈ {x i01 , . . . , x i04 } satisfy R |{xi 0 1,...,xi 0 4,x ′ } = {01101, 10011, 11000} or R |{xi 0 1,...,xi 0 4,x ′ } = {01100, 10010, 11001}. But this is a contradiction, since then R |{xi 0 1,...,xi 0 4,x ′ } is not closed under ternary majority, and so R / ∈ ID 2 . Thus the transformations end with R = R p 4 (X i0 ) (no other x j can be left). To conclude, from R we freezingly implemented R p 4 , and we are done. ✷ Theorem 19 Given any set of relations Γ such that Γ = ID 2 , then Γ f r equals one of the 13 frozen partial coclones in Definition 17.

It is easy to prove that the 13 frozen partial co-clones in Definition 17 are all distinct. Given two sets of relations Γ 1 and Γ 2 such that Γ 1 = Γ 2 = ID 2 , then to prove that Γ 1 f r = Γ 2 f r it is sufficient (according to Theorem 12) to show that f rP ol(Γ 2 ) = f rP ol(Γ 1 ).

The operations that we use to separate the different frozen partial co-clones are all ternary minority operations that are undefined on certain tuples. We denote these minority operations by mU (t 1 , . . . , t n ) where t 1 , . . . , t n are the tuples on which the minority operation is undefined. For example, Γ