Boundary control of partial differential equations using frequency domain optimization techniques

TL;DR

This paper introduces a frequency domain $H_$-control method for boundary control of PDEs, ensuring practical controller implementation and demonstrating effectiveness on reaction-diffusion and wave equations.

Contribution

The paper develops a novel frequency domain $H_$-control approach for boundary PDE control with practical, simple controllers, applicable to parabolic and hyperbolic systems.

Findings

Effective control of reaction-diffusion with input delay

Successful boundary anti-damping control of wave equation

Controllers are physically implementable and structurally simple

Abstract

We present a frequency domain based $H_\infty$-control strategy to solve boundary control problems for systems governed by parabolic or hyperbolic partial differential equation, where controllers are constrained to be physically implementable and of simple structure suited for practical applications. The efficiency of our technique is demonstrated by controlling a reaction-diffusion equation with input delay, and a wave equation with boundary anti-damping.

Figures

Equations

$$\Gamma :\{\begin{array}{ccc}\dot{x}& =& Ax\\ Px& =& u\\ y& =& Cx\end{array}$$

$${\Gamma }_s:\{\begin{array}{ccc}sx(s)& =& Ax(s)\\ Px(s)& =& u(s)\\ y(s)& =& Cx(s)\end{array}$$

$$\displaystyle\begin{array}[]{ll}\mbox{minimize}&\|T_{wz}(P,K)\|_{\infty}\\ \mbox{subject to}&K\mbox{ stabilizes }G\mbox{ in closed loop}\\ &K\in\mathscr{K}\end{array}$$

$$F(s)=\left[\begin{array}{cc}I& G(s)\\ -K(s)& I\end{array}\right],f(s)=\det F(s)$$

$$T=\left[\begin{array}{cc}{\left(I+KG\right)}^{-1}& -K{\left(I+GK\right)}^{-1}\\ {\left(I+GK\right)}^{-1}G& {\left(I+GK\right)}^{-1}\end{array}\right],$$

$$L[{\omega }_i,{\omega }_{i+1}]({\omega }_{i+1}-{\omega }_i)<|\tilde{f}(j{\omega }_i)|+|\tilde{f}(j{\omega }_{i+1})|.$$

$$L[{\omega }_i,{\omega }_{i+1}]({\omega }_{i+1}-{\omega }_i)<2{\gamma }^{\ast }+2\vartheta -\varphi ({\omega }_i)-\varphi ({\omega }_{i+1}),$$

$$e_1={\int }_{\overline{\omega }}^{\infty }f(\omega )d\omega \le \vartheta /2.$$

$$\frac{\overline{\omega }}{4}({\omega }_{i+1}-{\omega }_i)L[{\omega }_i,{\omega }_{i+1}]\le \vartheta /2,$$

$$e=\left|{\int }_0^{\infty }f(\omega )d\omega -{\int }_0^{\infty }P(\omega )d\omega \right|<\vartheta .$$

$$x_t(\xi ,t)=x_{\xi \xi }(\xi ,t)+c(\xi )x(\xi ,t),t\ge 0,\xi \in [0,L],$$

$$x(0,t)=0,$$

$$x(L,t)=u(t-D)$$

$$G(s,\xi )=\frac{x(\xi ,s)}{u(s)}=e^{-s}\frac{e^{\sqrt{s-\frac{1}{2}}\xi }-e^{-\sqrt{s-\frac{1}{2}}\xi }}{e^{\sqrt{s-\frac{1}{2}}L}-e^{-\sqrt{s-\frac{1}{2}}L}},$$

$$P_{\mathrm{red}}(s):\{\begin{array}{ccccccc}{\dot{x}}_r& =& A_r{x}_r& +& d& +& B_2{u}_r\\ z_r& =& C_1{x}_r& & & +& D_{12}{u}_r\\ y_r& =& C_2{x}_r& +& D_{21}d& & \end{array}$$

$$\displaystyle\begin{array}[]{lll}\mbox{minimize}&\|T_{d\to z_{r}}(P_{\rm red},K)\|_{2}\\ \mbox{subject to}&K(\mathbf{x})\mbox{ stabilizes }G_{\rm red}\mbox{ internally}\\ &K(\mathbf{x})\in\mathscr{K}_{2}\end{array}$$

$$K_0^{11}$$

$$K_0^{12}$$

$$K_0^{13}$$

$$K_0^{14}$$

$$K_0^{15}$$

$$\displaystyle\begin{array}[]{ll}\mbox{minimize}&\|(I+G_{\rm red}K_{0})^{-1}G_{\rm red}K_{0}-(I+GK)^{-1}GK\|_{\infty}\\ \mbox{subject to}&K\mbox{ stabilizes }G\mbox{ internally}\\ &K\in\mathscr{K}_{2}\end{array}$$

$$K_{11}^{\ast }$$

$$K_{12}^{\ast }$$

$$K_{13}^{\ast }$$

$$K_{14}^{\ast }$$

$$K_{15}^{\ast }$$

$$\displaystyle\begin{array}[]{ll}\mbox{minimize}&\left\|\left[\begin{array}[]{c}W_{e}(s)(I+G(s)K(s))^{-1}\\ W_{u}(s)K(s)(I+G(s)K(s))^{-1}\end{array}\right]\right\|_{\infty}\\ \mbox{subject to}&\mbox{$K$ stabilizes $G$ internally}\\ &K\in\mathscr{K}_{2}\end{array}$$

$$K_{11}^{\ast }$$

$$K_{12}^{\ast }$$

$$K_{13}^{\ast }$$

$$K_{14}^{\ast }$$

$$K_{15}^{\ast }$$

$$K_{11}^{\ast \ast }$$

$$K_{12}^{\ast \ast }$$

$$K_{13}^{\ast \ast }$$

$$K_{14}^{\ast \ast }$$

$$K_{15}^{\ast \ast }$$

$$x_{tt}(\xi ,t)$$

Full text

Keywords: Boundary control of PDEs, frequency-domain design, convection-diffusion, wave equation, $H_{\infty}$, structured feedback, infinite-dimensional systems

Pierre Apkarian is with ONERA, 2 Av. Ed. Belin, 31055, Toulouse, France Pierre.Apkarian@onera.fr

Dominikus Noll is with Université Paul Sabatier, Institut de Mathématiques, 118, route de Narbonne, 31062 Toulouse, France dominikus.noll@math.univ-toulouse.fr

1. Introduction

A recurrent issue in system control is whether, or to what extent, frequency-domain based $H_{\infty}$-control strategies originally developed for real-rational systems expand to infinite-dimensional processes. Success in rendering $H_{\infty}$-optimization fit to provide practically implementable controllers for infinite-dimensional systems should substantially foster the acceptance of PDE-modeling as a tool for control.

In response to this quest, we present a frequency domain based method to control infinite-dimensional LTI-systems, which is in particular suited for $H_{\infty}$-boundary control of parabolic and hyperbolic partial differential equations. Our method leads to practically implementable structured output feedback controllers for PDEs in such a way that the typical work-flow in control design is respected.

After briefly outlining our method, we will apply it in more detail to two infinite-dimensional $H_{\infty}$-control problems: boundary control of a reaction-diffusion system with input delay, as discussed in [1], and boundary control of an anti-stable wave equation to control noise and disturbance effects on duct combustion dynamics in a drilling pipe system [2, p.6], [3]. While the first study leads to a parabolic equation of retarded type, the second study leads to a system of neutral type, which poses new challenges to our frequency approach.

The structure of the paper is as follows. In section 2 we give the principal steps of our method. Stability is discussed in section 3, the role of the Nyquist test in optimization in section 3.1, its implementation in section 3.2. Sampling for performance is addressed in section 4. In section 6, we discuss the application of our method to a reaction-diffusion equation, and in section 7 to a wave equation.

2. Outline of the method

We start out with an infinite-dimensional LTI-system represented by a transfer function $G(s)$ with $p$ inputs and $m$ outputs, assumed well-posed in the sense of Salamon-Weiss [4, 5, 6, 7]. As principal application, we consider the case of a linearized parabolic or hyperbolic boundary control problem in state-space form

[TABLE]

with operators $A\in L(X,H),P\in L(X,\mathbb{R}^{p}),C\in L(X,\mathbb{R}^{m})$ on Hilbert spaces $H,X$ and finite-dimensional input and output spaces, where $X$ is densely embedded in $H$. Then under natural assumptions specified in [6, Sect. 2] the transfer function $G(s)$ of (1) is well-posed and obtained by applying the Laplace transform to (1), where every function evaluation $G(s)$ requires solving an elliptic boundary control problem

[TABLE]

Well-posedness means that $G(s)$ is holomorphic on a half-plane Re$(s)>\sigma$, but it may be convenient to require a little more, namely, that $G(s)$ extends meromorphically over a domain containing $\overline{\mathbb{C}}^{+}$. This is satisfied in all cases of practical interest, and guaranteed theoretically e.g. when $G$ is exponentially input/output stabilizable, see [8, Lemma 8.2.9 (i)(b), (ii)]. The meromorphic form of the transfer function is a necessary requirement for applicability of the Nyquist stability test.

After embedding $G(s)$ in a plant $P(s)$ with one or several closed-loop performance and robustness channels $T_{wz}(P,K)$, we set up the infinite-dimensional $H_{\infty}$-optimization problem

[TABLE]

where $\mathscr{K}$ represents a suitably chosen class of structured controllers with $m$ inputs and $p$ outputs. In this work we understand the term structured in the sense that controllers $K(\mathbf{x},s)$ depend differentiably on a vector $\mathbf{x}\in\mathbb{R}^{n}$ of tunable parameters, and have well-posed transfer functions $K_{ij}(\mathbf{x},s)$, typically with quasi-polynomial numerators and denominators. Such control laws combine real-rational elements with input and output delays, and can therefore be physically implemented. In the optimization procedure it will also be necessary to know the finite number of rhp poles of $K(\mathbf{x})$ for every $\mathbf{x}$.

Since (6) as a rule cannot be solved exactly, we use an inexact bundle trust-region method as in [9, 10, 11], which guarantees stability of the closed loop, and approximates $H_{\infty}$-performance up to a user specified precision. The following scheme presents our method in a more formal way.

steps of this scheme require further explanations, which we provide in the following sections.

3. Stability test

Let us recall that with the definitions

[TABLE]

the inverse $T(s)=F(s)^{-1}$ is given as

[TABLE]

and we call the closed-loop system $(G,K)$ stable in the $H_{\infty}$-sense, or simply stable, if the transfer function $T$ belong to the Hardy space ${\bf H}_{\infty}(\mathbb{C}^{+},\mathbb{C}^{(m+p)\times(m+p)})$. As is well-known, $H_{\infty}$-stability is equivalent to the absence of unstable poles in tandem with boundedness of $T(s)$ on $j\mathbb{R}$. We are interested in situations, where absence of unstable poles of $T(s)$ can be verified by the Nyquist stability test.

Systems arising from parabolic equations are of retarded type and typically satisfy the spectrum decomposition assumption, which means that they have only a finite number of unstable poles. The Nyquist stability test may therefore be applied directly, to the effect that in order to guarantee absence of unstable poles in the loop $f(j\omega)=\det(I+G(j\omega)K(\mathbf{x}^{0},j\omega))$ has to wind $n_{p}$ times around the origin in the clockwise sense, where $n_{p}$ the number of rhp poles of $G$ and $K(\mathbf{x})$ together.

In order to address the case where $f$ has a finite number of poles on $j\mathbb{R}$, we consider the following construction, which avoids the usual $\epsilon$-indentation of the Nyquist contour into the rhp. We choose a holomorphic function $h$ on a domain containing $\overline{\mathbb{C}}^{+}$ such that $h(s)\not=0$ on $\mathbb{C}^{+}$, $\lim_{s\to\infty}h(s)=1$ on $\overline{\mathbb{C}}^{+}$, and such that $h$ has a zero of order $p$ at $\pm j\omega$ precisely when $F$ has a pole of order $p$ at $\pm j\omega$. Let $\widetilde{f}=fh$, $\operatorname{\mathsf{D}}$ a Nyquist D-contour into the rhp with $[-j\overline{\omega},j\overline{\omega}]\subset\operatorname{\mathsf{D}}$ containing in its interior all rhp poles of $F$. Then the modified Nyquist curve $\widetilde{f}\circ\operatorname{\mathsf{D}}$ has the same winding number as the original Nyquist curve $f\circ\operatorname{\mathsf{D}}_{\epsilon}$ with sufficiently small $\epsilon$-indentations.

From the moment onward a controller $K(\mathbf{x})$ has been identified closed-loop stabilizing using the Nyquist test, the nonsmooth optimization method, when considering a trial step $K(\mathbf{x}+d\mathbf{x})$ away from the current iterate $\mathbf{x}$, will re-compute the winding number to check stability of the the loop with $K(\mathbf{x}+d\mathbf{x})$. In those cases where the number of poles of $K(\mathbf{x})$ is independent of $\mathbf{x}$, this means we simply have to assure that the winding number $n_{p}$ does not change as we go from $\mathbf{x}$ to $\mathbf{x}+d\mathbf{x}$, which requires preventing the Nyquist curve from crossing the origin. Stability in the $H_{\infty}$-sense for the new $K(\mathbf{x}+d\mathbf{x})$ will then follow under the proviso that the closed-loop transfer function $T(\mathbf{x}+d\mathbf{x})$ remains bounded on $j\mathbb{R}$, which is the case when $G,K(\mathbf{x})$ are bounded on $j\mathbb{R}$, and occurs in particular if these transfer functions are proper. With these preparations the situation for parabolic systems is covered by the following:

Theorem 1**.**

Suppose process $G(s)$ and controller $K(s)$ are well-posed, extend meromorphically into a domain containing $\overline{\mathbb{C}}^{+}$, and satisfy the following conditions:

1. i.

$F$ has no zeros on $j\mathbb{R}$ and only a finite number $N_{p}$ of poles in $\overline{\mathbb{C}}^{+}$, $n_{p}$ of which are in $\mathbb{C}^{+}$. 2. ii.

There are no pole/zero cancellations on $j\mathbb{R}$. 3. iii.

There exist a frequency $\overline{\omega}>0$ and $\alpha>0$ such that Re$f(j\omega)>\alpha$ for all $\omega\in[\overline{\omega},\infty)$. 4. iv.

$G,K$ are bounded on $j\mathbb{R}\setminus[-j\overline{\omega},j\overline{\omega}]$. 5. v.

$G,K$ have strongly (exponentially) stabilizable and detectable state-space realizations.

Suppose the modified Nyquist curve $\widetilde{f}\circ\operatorname{\mathsf{D}}$ winds $n_{p}$ times around the origin in the clockwise sense. Then the closed-loop $T(s)$ is strongly (exponentially) stable.

• $\square$*

The theorem was proved in [12] under the more standard assumptions that $G,K$ are proper and $\lim_{s\to\infty}f(s)\not=0$ on $\overline{\mathbb{C}}^{+}$ exists, but since our present statement concerns only a finite Nyquist contour, $\widetilde{f}\circ\operatorname{\mathsf{D}}$, the proof can be adapted with minor changes.

The full force of Theorem 1 is needed when it comes to dealing with hyperbolic systems. Here the situation is complicated because in open loop neutral systems may have infinitely many rhp poles in a strip $0\leq{\rm Re}(s)<\alpha$, so condition (iii) may fail, even though $G$ may still be well-posed. When this is the case, it is impossible to use the Nyquist test directly even if an initial stabilizing controller $K_{0}=K(\mathbf{x}^{0})$ is given, because the Nyquist curve $f(s)=\det(I+G(s)K_{0}(s))$ winds infinitely often around the origin. In this event the method explained in the following section is helpful.

3.1. Enabling the Nyquist test

Assuming that is has been verified by some other means that the closed-loop system $G_{0}=G(I+K(x_{0})G)^{-1}=:{\tt feedback}(G,K(\mathbf{x}^{0}))$ is indeed stable. Then a small gain argument tells us that ${\tt feedback}(G_{0},K)$ remains stable for stable controllers $K$ satisfying $\|K\|_{\infty}<1/\|G_{0}\|_{\infty}$. Since ${\tt feedback}(G_{0},K)=$ ${\tt feedback}({\tt feedback}(G,K(\mathbf{x}^{0})),K)$ $={\tt feedback}(G,K(\mathbf{x}^{0})+K)$, we see by letting $K$ the difference $K(\mathbf{x})-K(\mathbf{x}^{0})$ that ${\tt feedback}(G,K(\mathbf{x}))$ is stable for $\|K(\mathbf{x})-K(\mathbf{x}^{0})\|_{\infty}<1/\|G_{0}\|$, and this can now be verified by applying the Nyquist test to $f_{0}(s)=\det(I+G_{0}(s)K(s))$, where $K=K(\mathbf{x})-K(\mathbf{x}^{0})$.

Here Theorem 1 applies indeed to $G_{0}$, $K$, because $G_{0}$ has no poles on $\overline{\mathbb{C}}^{+}$, so that i.-iv. are satisfied for $G_{0},K$ provided iii. was from start satisfied for $G,K_{0}$. We have, however, to recall that despite stability of the loop, $f_{0}$ will typically not have a limit as $s\to\infty$ on $\overline{\mathbb{C}}^{+}$, so we will rely on condition iv., which assures that outside the band $[-\overline{\omega},\overline{\omega}]$ the Nyquist curve is in no danger to turning around 0.

From here on, we can proceed just as in the previous case for retarded systems, where now $N_{p}=n_{p}$ is the known number of rhp poles of $K=K(\mathbf{x})-K(\mathbf{x}^{0})$.

Remark 1**.**

Should our initial stabilizing controller $K_{0}$ be unstable, it may be preferable to use controllers of the form $K_{0}+K(\mathbf{x})$ for $K(\mathbf{x})$ stable, which gives rise to a modified structure.

The question remains in what sense a state-space representation of $T$ given by (7) will be stable. This is decided by the following

Theorem 2**.**

Suppose $G,K$ are well-posed transfer functions which admit strongly (exponentially) stabilizable state-space realizations. Suppose the closed loop transfer function satisfies $T\in{\bf H}_{\infty}$. Then the generator of the state-space representation of the closed loop is strongly (exponentially) stable.

Proof.

For exponential stability, according to Morris [13, Theorem 5.2] it suffices to show that the closed loop is exponentially stabilizable and exponentially detectable. By Staffans [8, Lemma 8.2.7] this follows as soon as each of the components $G,K$ is individually exponentially stabilizable and detectable, but the latter is true by hypothesis.

For the statement concerning strong stability we use [8, Lemma 8.2.7] again, which now guarantees that the closed loop $T(G,K)$ is strongly stabilizable and detectable. Since it is $H_{\infty}$-stable, we can invoke [8, Theorem 8.2.11 (ii)] to infer that the closed loop is also strongly stable. ∎

3.2. Sampling for the Nyquist test

During optimization, the Nyquist test is applied at every candidate controller $K(\mathbf{x})$ to check stability of the loop ${\tt feedback}(G,K(\mathbf{x}))$, respectively for neutral systems, of the preliminary stabilized system ${\tt feedback}(G_{0},K(\mathbf{x})-K(\mathbf{x}^{0}))$. For that we have to sample the modified Nyquist curve $\widetilde{f}=fh$ at frequencies $0=\omega_{0}<\dots<\omega_{N}=\overline{\omega}$ such that the winding number of $\widetilde{f}\circ\operatorname{\mathsf{D}}$, based on a typical D-contour extending into the rhp, coincides with the winding number of the closed polygon $P_{\widetilde{f}}$ obtained as $\widetilde{f}(j\omega_{0}),\dots,\widetilde{f}(j\omega_{N})$, $\widetilde{f}(-j\omega_{N})$, $\widetilde{f}(-j\omega_{N-1}),\dots,\widetilde{f}(-j\omega_{0})$. The latter can then be computed conveniently using the ray crossing algorithm (see Fig. 1)

To begin with, we need a cutoff frequency $\overline{\omega}>0$ such that every Nyquist contour $\operatorname{\mathsf{D}}$ containing the segment $[-j\overline{\omega},j\overline{\omega}]$ contains all $n_{p}$ unstable poles of $G,K$ in its interior. We then sample on the segment $[0,\overline{\omega}]$, and replace the curved part of the D-contour $\widetilde{f}\,\circ$ 

$\operatorname{\mathsf{C}}$

by the segment $[\widetilde{f}(-j\overline{\omega}),\widetilde{f}(j\overline{\omega})]$. For the latter to be authorized, the closed curve $\gamma$ concatenated by this segment and $\widetilde{f}\,\circ$ 

$\operatorname{\mathsf{C}}$

has to satisfy ${\rm ind}(\gamma,0)=0$ and should not contain any of the $n_{p}$ poles of $F$ in its interior. This is for instance the case if we choose $\overline{\omega}$ as in condition iv. of Theorem 1. Then we sample $\widetilde{f}$ on $[0,\overline{\omega}]$ such that none of the curves $\gamma_{i}$ obtained by concatenating the segment $[\widetilde{f}(j\omega_{i+1}),\widetilde{f}(j\omega_{i})]$ with the piece $\widetilde{f}([j\omega_{i},j\omega_{i+1}])$ of the modified Nyquist curve encircles [math]. This can be assured my the method in [12]. Call a mapping $L[\cdot,\cdot]:\mathbb{R}^{2}\to\mathbb{R}$ a first-order bound of $\widetilde{f}$ if $|\widetilde{f}^{\prime}(j\omega)|\leq L[\omega^{-},\omega^{+}]$ for all $\omega\in[\omega^{-},\omega^{+}]$. Then we have the following

Lemma 1**.**

Suppose that for fixed $K$ the cutoff frequency $\overline{\omega}>0$ is as above, and the sampling nodes satisfy the condition

[TABLE]

Then the winding numbers of the modified Nyquist curve $\widetilde{f}\circ\operatorname{\mathsf{D}}$ and the approximating closed polygon $P_{\widetilde{f}}$ agree. **

For the proof see [12, sect 3.]. $\square$

4. Sampling for performance

Sampling for $H_{\infty}$-performance was also analyzed in [12] and can again be based on a first-order bound $L[\cdot,\cdot]$, now for the function $\phi(\omega)=\overline{\sigma}\left(T_{wz}(P(j\omega),K(j\omega))\right)$. We recall the following

Lemma 2**.**

Let $\gamma^{*}=\max\{\phi(\omega_{i}):i=1,\dots,N\}$ for a given controller $K$ and a corresponding sampling $\omega_{i}$. Let $\vartheta>0$ be a user-specified tolerance. If the nodes $\omega_{i}$ satisfy

[TABLE]

then the true $H_{\infty}$ norm is within tolerance $\vartheta$ of its estimated value, that is, $\gamma^{*}\leq\|T_{wz}(P,K)\|_{\infty}\leq\gamma^{*}+\vartheta$. $\square$ **

For the proof see [12, sect. 5]. In that work we have compared sampling for stability via (8) and sampling for $H_{\infty}$-performance based on (9) on a large test bench including finite and infinite dimensional systems. The results fairly consistently show that performance requires at least 10 times more nodes $\omega_{i}$ than sampling to assure that the Nyquist test is correct. This leads to the following significative meta-theorem: $H_{\infty}$-performance is 10 times more costly than mere stability.

In those cases where a channel $T=T_{wz}(P,K)$ for $H_{2}$-optimization is available, we need a sampling $\omega_{i}$ for $f(\omega)={\rm trace}\left(T(j\omega)^{H}T(j\omega)\right)$. Assume that a first-order bound $L[\omega^{-},\omega^{+}]$ for $f$ is available, and let $P$ be the piecewise linear function corresponding to the polygon with nodes $(\omega_{i},f(\omega_{i}))$, $0=\omega_{0}<\dots<\omega_{N}=\overline{\omega}$ with $P(\omega)=0$ for $\omega>\overline{\omega}$.

Lemma 3**.**

Let $\vartheta>0$ be a user specified tolerance and suppose the cutoff frequency $\overline{\omega}>0$ is such that

[TABLE]

Suppose the interval $[0,\overline{\omega}]$ is sampled with nodes $\omega_{i}$ such that

[TABLE]

then the error satisfies

[TABLE]

Proof.

Let $e_{1}$ be the error of the high frequency contribution satisfying $e_{1}<\vartheta/2$. Now the error of the low frequency part is $e_{1}=\left|\int_{0}^{\overline{\omega}}f(\omega)d\omega-\int_{0}^{\overline{\omega}}P(\omega)d\omega\right|\leq\sum_{i=0}^{N-1}\frac{1}{4}(\omega_{i+1}-\omega_{i})^{2}L[\omega_{i},\omega_{i+1}]\leq\sum_{i=0}^{N-1}(\omega_{i+1}-\omega_{i})\vartheta/2\overline{\omega}=\vartheta/2.$ Hence altogether $e=e_{1}+e_{2}<\vartheta$. ∎

Remark 2**.**

It is clear that (10) is much more binding than (9), because sampling to assure the exactness of the maximum value within a tolerance $\vartheta>0$ has only to be precise at frequencies close to the maximum, whereas approximating an integral requires good approximation on the whole $[0,\overline{\omega}]$. This suggests avoiding $H_{2}$-optimization if possible. Since robustness requirements further press to avoid $H_{2}$-optimization, we presently seek for workarounds.

5. Optimization

With a computable test for stability and a method to approximate the objective function $\|T_{wz}\|_{\infty}$ available, we run our nonsmooth optimization method based on [12, 10], with the interpretation of inexact function and subgradient evaluations as in [9]. For this we have to recall subgradient evaluation as discussed in [14]. Suppose $\omega_{i}$ is one of the sample frequencies where the maximum $\gamma^{*}=\phi(\omega_{i})$ of the approximation is attained with error $\gamma^{*}\leq\|T_{wz}(\mathbf{x})\|_{\infty}\leq\gamma^{*}+\vartheta$. Then an approximate subgradient is generated by computing one or several maximum eigenvectors $\Phi_{i}$ of $T_{wz}(\mathbf{x},\omega_{i})^{H}T_{wz}(\mathbf{x},\omega_{i})$ and using formulas (12) or (13) of [14, Sect. IV]. Inspecting those shows that for the nearby frequency $\omega$ where $\|T_{wz}(\mathbf{x})\|_{\infty}$ is in reality attained, the mismatch between the estimated eigenvector $\Phi_{i}$ and one of the true maximum eigenvectors $\Phi$ of the $H_{\infty}$-norm at $\omega$ is proportional to the error in the function values, but with proportionality constant depending on the reciprocal of the eigenvalue gap at $T_{wz}(\mathbf{x},\omega)$. Since only finitely many frequencies are active, the eigenvalue gap cannot become arbitrarily small. From the same reason it remains bounded in the neighborhood of any of the accumulation points of the sequence of serious iterates generated by the bundle or the bundle trust-region method. That suggests indeed an interpretation of our method as an instance of the inexact bundle trust-region method.

For a recent thorough convergence analysis of the bundle method with inexact function and subgradient evaluations in an infinite-dimensional setting we refer to Hertlein and Ulbrich [15].

Starting the algorithm as presented in [12, 10] at the closed-loop stabilizing $K_{0}=K(\mathbf{x}^{0})$, the Nyquist test is used at every new iterate $K(\mathbf{x}^{+})$ to check whether the loop ${\tt feedback}(G,K(\mathbf{x}^{+}))$, respectively ${\tt feedback}(G_{0},K(\mathbf{x}^{+})-K(\mathbf{x}))$ for the case of neutral systems, is stable. If this is not the case, a backtracking step $\mathbf{x}_{\alpha}=\mathbf{x}+\alpha(\mathbf{x}^{+}-\mathbf{x})$ for $0<\alpha<1$ is made such that $K(\mathbf{x}_{\alpha})$ is still stabilizing, and a repelling cutting plane is included in the bundle, using e.g. the closed-loop sensitivity function $S(\mathbf{x})=\|(I+GK(\mathbf{x}))^{-1}\|_{\infty}$ as a stability barrier function. See [12, sect. 4] for details.

Remark 3**.**

A special situation occurs if $G,K$ are stable, or if $G_{0}$ instead of $G$ is used as described in section 3.1. Here the Nyquist curve $f(j\omega)$ does not at all wind around the origin, and as a rule stays outside a conical or parabolic region $R_{\alpha,r}=\{s\in\mathbb{C}:{\rm Im}(s)^{2}<\alpha(r-{\rm Re}(s))\}$ for certain $\alpha,r>0$. The constraint ${\rm Im}(f(\mathbf{x},j\omega))^{2}\geq\alpha(r-{\rm Re}(f(\mathbf{x},j\omega)))$ can be easily included in the constrained program (6). This is more reliable than just preventing $f(\mathbf{x},j\omega)$ from crossing the origin. Note that even when the Nyquist test can be applied directly to $G,K(\mathbf{x})$, e.g. for parabolic or first-order hyperbolic systems, it may for numerical reasons be interesting to apply it in just the same way to $G_{0},K$.

A more conservative but robust way to address the same problem is to add a disk-margin constraint of the form $\|S\|_{\infty}\leq 1/\alpha$.

6. Output-feedback control of a reaction-diffusion equation with input delay

In this section we discuss a one-dimensional reaction diffusion equation with delayed Dirichlet boundary control

[TABLE]

where $x(\cdot,t)$ denotes the state of the system, $u(t)$ the control, $D$ the delay, and where we assume that a finite number of measured outputs $y_{1}(t)=x(\xi_{1},t),\dots,y_{m}(t)=x(\xi_{m},t)$ at sensor positions $\xi_{i}\in[0,L]$ are available for control. A similar control scenario is discussed in Prieur and Trélat [1] under the assumption of full state measurement. Related work is for instance Sano [16], where $H_{\infty}$-control of a heat exchanger is discussed, or [10], where a reaction-convection-diffusion equation with simultaneous boundary and distributed control and a van de Vusse reactor of a coupled system of reaction convection-diffusion equations again with combined boundary and distributed control are discussed without input delay, but with a single point measurement as output.

In the present study, we strive to control the system with a finite-dimensional output feedback controller $K(\mathbf{x})$ of simple structure, which could conveniently be implemented, and yet gives satisfactory performances in closed loop.

Performance specifications of the reaction-diffusion PDE are chosen so that responses to non-zero initial conditions show reasonable behavior in terms of damping and settling time. This could be addressed by an $H_{2}$-performance specification, but as we show may also be successfully controlled by way of suitably chosen $H_{\infty}$-specification. The latter is advantageous as soon as additional robustness aspects of the design are called for.

Working incrementally, and starting with the case of a single measurement at the mid-point $\xi=L/2$, our analysis indicates that $5$ equidistant measurements are enough to achieve good responses against initial conditions, while mere stability could be assured even on the basis of a single measurement e.g. at $\xi=0$.

In our numerical testing we adopt the choices $L=2\pi$, $D=1$ and $c(x)=\frac{1}{2}$ from [1], where the open-loop $F(s)$ has one unstable pole at $s=\frac{1}{2}$, and an infinity of stable double poles at $s_{k}=\frac{1}{2}-\frac{k^{2}\pi^{2}}{L^{2}}$ following a retarded pattern. This is understood, as the semi-group of the equation is sectorial [8, p. 150]. As indicated in previous chapters, this allows direct application of the Nyquist test to check $H_{\infty}$-stability of the closed loop.

Remark 4**.**

Well-posedness of the system (6) in the sense of Salamon-Weiss can be deduced from the functional analytic setting in [1], or from the general approach in [6].

The transfer function of (6) can be computed analytically as

[TABLE]

and for $L=2\pi$ the system has one unstable pole.

Proposition 1**.**

Suppose a finite-dimensional structured controller $K(\mathbf{x})$ with $m$ inputs and $p=1$ output is found which stabilizes system* (6) internally in the $H_{\infty}$-sense. Then the closed loop is even exponentially stable.*

Proof.

By theorem 2 this follows as soon as each of the components $G,K$ is individually exponentially stabilizable and detectable. Since $K$ is finite-dimensional it has clearly an exponentially stabilizable and detectable state-space model. For $G$ exponential stabilizability may be deduced from [1], because the infinite-dimensional state-feedback controller the authors construct has the same control input as (6). Exponential detectability on the other hand follows from the fact that the differential operator in (6) is self adjoint, and that in the adjoint system the five outputs are turned into 5 inputs, one of which is the same as the single input in (6), but now without the delay. Exponential detectability therefore follows from the fact that (6) is exponentially stabilizable without the delay. ∎

This still leaves the problem of finding a preliminary stabilizing controller $K_{0}=K(\mathbf{x}^{0})$ of the pre-defined structure. As we indicated in previous sections, the latter as a rule requires heuristic methods even for very simple structures. The advantage we have in the case of the present parabolic study (6) is that we can check via the Nyquist test whether a given controller is stabilizing.

6.1. Model matching approach

Model matching is a sophisticated control scenario, where specifications are pursued indirectly. It is covered by algorithm 1, and we believe it is particularly suited for PDE-control, where models of different grid scales arise naturally. Here we use model matching to address the reaction of the system to a non-zero initial value. The method consists in two steps (a) and (b).

Non-vanishing initial values may be regarded as disturbances $d$ acting on the system state as in Fig. 2. We now assume that we have to regulate against functions $x_{0}(\xi)\not=0$ of a given bandwidth. In a first step (a) we therefore compute a reduced state-space model $G_{\rm red}(s)$ of $G(s)$, in which the resolution of the state $x_{r}(\xi)$ reflects the resolution of the potential $x_{0}(\xi)$ accurately. In the present study we regulate against initial conditions with resolution comparable to that considered in [1], which leads us to a finite-difference discretization of (6) with $50$ spatial steps or states, complemented by a $3$rd-order Padé approximation of the input delay adding $3$ more states. This coarse grid model $G_{\rm red}$ is embedded into a plant $P_{\rm red}$ expressing control requirements in terms of damping and settling time in responses to initial conditions. Here this consists in optimizing the root mean-square energy value of the output signal $z_{r}=(z_{1},z_{2})$ in response to the white noise disturbance $d$ on the state $x_{r}$:

[TABLE]

where $x_{r}\in\mathbb{R}^{54}$ is the reduced state, $d$ is the exogenous input, understood to represent the impulse caused by the non-zero initial value, and $z_{r}=(W_{x}x_{r},W_{u}u)$ is regulated similar to what is used in LQG-control, with filters $W_{x}=I$ and $W_{u}(s)=\frac{s}{1+s/a}$, $a=100$, the latter adding another state to $x_{r}$. The reduced output $y_{r}\in\mathbb{R}^{5}$ of $G_{\rm red}$ represents the 5 distributed measurements (see Fig. 5) in the coarse finite-difference discretization. Our testing shows that five equidistant measurements along $[0,\,L]$ are sufficient to achieve well behaved responses to initial conditions.

We then solve the structured $H_{2}$-optimization problem for the channel $d\to z_{r}$ of $P_{\rm red}$ in Fig. 2, using systune [17] based on [14, 18]. Trial and error reveals the interesting fact that fairly low-order controllers $K\in\mathscr{K}_{2}$ of order 2 are adequate, leading to $\mathbf{x}\in\mathbb{R}^{21}$ for $K$:

[TABLE]

where the closed-loop transfer $T_{d\to z_{r}}(P_{\rm red},K)$ is the lower linear fractional transformation $F_{l}(P_{\rm red},K)$. The solution of (15) is $K_{0}$ in Figure 3 right. The difference with traditional $H_{2}$-control is that $K$ does not have observer structure, but the structure $K\in\mathscr{K}_{2}$ we imposed. The resulting controller $K_{0}\in\mathscr{K}_{2}$ is obtained as

[TABLE]

A simulation for initial condition $x_{0}(\xi)=\xi(L-\xi)$ is shown in Fig. 5 (left). While this produces the expected good results for $G_{\rm red}$, we now have to check whether $K_{0}$ also stabilizes the infinite-dimensional system $G(s)$. It turns out that this is the case, as the Nyquist test reveals, so that we now proceed to the second part (b) of the model matching method, where the controller is further optimized with regard to the full model.

As a result of step (a) of the model matching procedure we have so far obtained a reference model, ${\tt feedback}(G_{\rm red},K_{0})$ and a controller $K_{0}$ of the desired structure which stabilizes $G_{\rm red}$. Application of the Nyquist test, in tandem with boundedness of the closed loop transfer function on $j\mathbb{R}$, show that $K_{0}$ also stabilizes the infinite dimensional system $G(s)$. Fig. 4 (left) shows one clockwise encirclement, computed by the ray-crossing algorithm, which due to the known single unstable pole in $G$, and absence of unstable poles in $K_{0}$, confirms the absence of unstable closed-loop poles. Taking into account that $G,K_{0}$ are proper shows that the closed loop transfer function is bounded on $j\mathbb{R}$, hence the loop is $H_{\infty}$-stable, and by Proposition 1, is exponentially stable.

Remark 5**.**

The fact that $K_{0}$ stabilizes not only $G_{\rm red}$, but even $G$, could be called accidental. However, recall that within most structures $K\in\mathscr{K}$, practical methods leading to a stable closed loop are necessarily heuristics, so remain equally accidental. What can be said in favor of our method to obtain $K_{0}$ is that it is the result of a local optimization procedure.

In a second step (b) of the model matching procedure, corresponding to step 3 of algorithm 1, the preliminary stabilizing controller $K_{0}$ is refined through an $H_{\infty}$ model matching problem shown in Fig. 3, which takes the true infinite-dimensional dynamics in $G(s)$ accurately into account. In this step, we have to solve an infinite-dimensional structured $H_{\infty}$-control problem covered by the general form (6),

[TABLE]

for which we use the bundle algorithm of [10, 9, 12], initialized at $K=K_{0}$. As a result of optimization the $H_{\infty}$-norm of the mismatch channel $r\to z=y-y_{r}$ is reduced from $1.81$ at $K_{0}$ to $0.84$ at $K^{*}$.

The optimized controller $K^{*}$ obtained is

[TABLE]

its simulation is shown in Figs. 5 (right) and the corresponding Nyquist plot in Fig. 4 (right). The results resemble those obtained by studies based on full-state information, see [1].

6.2. Mixed sensitivity approach

A classical approach to improve system performance is via minimization of the sensitivity function $S:=(I+GK)^{-1}$. Using appropriate weighting filters $W_{e}(s)$, one can achieve better transient and steady-state responses for all references signals of finite energy. Introducing a penalization of the control effort as before leads to a mixed-sensitivity design problem:

[TABLE]

which is a particular case of (6).

Actuation takes place at the edge $\xi=L$, and its effect propagates with a unit delay $e^{-s}$ along the spatial dimension. The settling time is essentially determined by the slow dynamics that correspond to the states farthest away from the actuation point. On the other end, states closer to $\xi=L$ settle faster, but in turn are affected by much more turbulent transients, as for instance seen in Fig. 6. This suggests shaping response surfaces using weightings which take the distance to $\xi=L$ into account. This leads us to $W_{e}={\rm diag}(c_{1},c_{2},c_{3},c_{4},c_{5})$ with individually adapted $c_{i}$, where for simplicity static weightings are sought. Note that $c_{1}$ corresponds to the edge $\xi=0$, while $c_{5}$ is associated with $\xi=L$. As before, penalization of the control effort uses a high-pass filter $W_{u}(s)=(s/10)/(1+s/1e3)$.

Fig. 6 (left) shows the simulation of $G$ in closed loop with a first optimal controller $K^{*}\in\mathscr{K}_{2}$ given in (6.2) obtained via nonsmooth optimization (27) started at $K_{0}$ from (32) and using $W_{e}^{(1)}={\rm diag}(3,0,0,0,0)$. States close to the edge $\xi=0$ have excellent settling times, but transient wobbles manifest themselves at the opposite end $\xi=L=2\pi$ (see Fig. 6 left).

[TABLE]

Increasing the cost at $\xi=L$ via $W_{e}^{(2)}={\rm diag}(1,0,0,0,0.2)$, and starting optimization at $K^{*}$, now leads to the optimal controller $K^{**}$ given in (6.2), which removes this undesirable effect. Simulation of the loop with $K^{**}$ including the control signal is shown in Fig. 6 (right), the Nyquist plot in 6 (right). The optimal controller is obtained as

[TABLE]

Finally, note that it is possible to obtain even faster responses by accepting more aggressive control signals and therefore obtaining a more academic than practical solution. By proposition 1 all controllers obtained are exponentially stabilizing.

Altogether this study shows that by way of program (6) it is possible to conveniently control the reaction-diffusion equation (6) with a single input with delay and 5 distributed measurements by synthesizing a finite-dimensional low-order controller such that the result matches the result obtained in [1] using full state feedback.

7. Control of an anti-stable wave equation

In this second study, we discuss the following boundary control system

[TABLE]

where $q>0$, $q\not=1$. The state of the system is $x(\cdot,t),x_{t}(\cdot,t)$, the control applied at the boundary $\xi=1$ is $u(t)$, and we assume that the measured outputs are

[TABLE]

The system has been discussed previously in [2], and [19, 3], where potential applications are mentioned. Its well-posedness can be seen from the functional analytic set-up in [19, 3], and from the general approach to well-posedness of 1D hyperbolic systems in [20].

The transfer function of (7) is obtained from the elliptic boundary value problems

[TABLE]

which in this particular situation can be solved analytically:

[TABLE]

From this general formula the transfer function of (7) is $G(s)=\left[G(0,s);G(1,s);sG(1,s)\right]=:\left[G_{1}(s);G_{2}(s);G_{3}(s)\right]$. The main challenge in the hyperbolic system (7) is that along with the unstable pole at $s=0$ it exhibits an infinite number of unstable poles on a line Re$(s)=\sigma>0$. This means that the Nyquist test is not directly applicable.

7.1. Preliminary stabilization

Following our scheme in algorithm 1, the first step is to provide a preliminary stabilizing controller $K_{0}=K(\mathbf{x}^{0})$ of a simple pre-defined structure. We have to stabilize the system

[TABLE]

where $Q=(1+q)/(1-q)$.

A first question is whether $G$ can be stabilized by a finite-dimensional controller. Ignoring the input $y_{2}$, which for stabilization is not required, we choose the structure $K(\mathbf{x})=\left[n_{1}(s)/d(s),0,n_{3}(s)/d(s)\right]$, with $\mathbf{x}$ gathering the unknown coefficients of the polynomials $n_{i}(s),d(s)$ with deg$(n_{i})\leq{\rm deg}(d)$. Stability of the closed loop $T(s)=G(s)/(1+G_{1}(s)K_{1}(s)+G_{3}(s)K_{3}(s)$ leads to testing whether the quasi-polynomial

[TABLE]

arising in the denominator of $T(s)$ is stable, i.e., has its roots in $\mathbb{C}^{-}$. While there exist general methods to check stability of quasi-polynomials, cf. [21], an ad hoc solution is here to choose $n_{3}=d$, whence the quasi-polynomial simplifies to

[TABLE]

where $c(s)=n_{1}(s)/(1-q)$ and deg$(n_{1})={\rm deg}(c)\leq{\rm deg}(d)$. If we choose $d(s)=s+x_{1}$ and $c(s)=x_{2}s+x_{3}$, then stability of the loop is equivalent to stability of the quasi-polynomial

[TABLE]

which is covered by the discussion in [22]. In their terminology we have $a_{0}=0$, $a_{1}=x_{1}$, $b_{2}=0$, $b_{1}=1$, $b_{0}=1$. We are then necessarily in the case $m=1$, $\mu_{0}=0$ of [22], so the quasi-polynomial $P(s)$ can only be stable if $x_{1}>-1$. Moreover the family $P_{h}(s)=A(s)+B(s)e^{-hs}$ is stable for all $0\leq h<h_{\sigma,0}$, where $h_{\sigma,0}>0$ is determined as follows. Let $\omega_{\sigma}$ be the positive real solution of

[TABLE]

and let $h_{\sigma,0}$ be the smallest positive solution $h$ of

[TABLE]

where $\arg(\cdot)\in[0,2\pi)$. If we let $x_{1}=1>-1$, then $4\omega_{\sigma}^{3}=1$, $\omega_{\sigma}=4^{-3}$, and we get

[TABLE]

and since our delay is $h=1$, this must now be solved for $x_{2},x_{3}$ so that $h_{\sigma,0}>1$. For instance $x_{2}=-1$ and $x_{3}=-4^{-3}$ gives argument $\pi/4$ in the formula, so that $h_{\sigma,0}=16\pi>1$. The leads to the finite-dimensional stabilizing controller $K_{0}$ for $G$:

[TABLE]

A second way to seek preliminary stabilization of (7) is to stick to the form $K=\left[\frac{n(s)}{d(s)}\;\;0\;\;1\right]$, but allow $n(s),d(s)$ to be quasi-polynomials, trying to simplify the denominator quasi-polynomial $P$ as much as possible. A very straightforward way is to let $d(s)=a(s)+e^{-s}b(s)$ with $a(s),b(s)$ polynomials, then the denominator quasi-polynomial simplifies to

[TABLE]

where $c(s)=n_{1}(s)/(1-q)$. If we now let $c(s)=-sb(s)$, then $K$ will be stabilizing in the $H_{\infty}$-sense as soon as $a(s)$ is stable, because the factor $s$ cancels with the factor $s$ in the numerator. If we choose $a(s)=s+c_{0}$, $b(s)=-c_{0}$ for some constant $c_{0}>0$, then we obtain the controller

[TABLE]

which in [3] was obtained using the back-stepping technique. Since only input and output delays along with real-rational terms arise, such controllers are implementable, so we are still in line with our general purpose of computing practically useful controllers.

7.2. Performance optimization

Let us now discuss a more systematic way which not only leads to preliminary stabilizing $K_{0}$, but also allows performance optimization. In order to compare with [3], we optimize again against the effect of non-zero initial values, using the output $y\in\mathbb{R}^{3}$, and aiming as before at a convenient implementable controller structure.

We start by putting the system $G$ in feedback with the controller $K_{0}=[0\;0\;1]$, which leads to $\widehat{G}=G/(1+G_{3})$, where

[TABLE]

Re-write this as $\widehat{G}=\widetilde{G}+\Phi$, where $\widetilde{G}$ is now real-rational and still unstable, while $\Phi$ gathers the infinite dimensional part, but is stable. Then we use that stability of the closed loop $(\widetilde{G}+\Phi,K)$ is equivalent to stability of the loop $(\widetilde{G},{\tt feedback}(K,\Phi))$, as explained in Fig. 8:

Then we construct a finite-dimensional structured controller $\widetilde{K}=\widetilde{K}(\mathbf{x})$ which stabilizes $\widetilde{G}$. The controller $K$ is then recovered from $\widetilde{K}$ through the equation $\widetilde{K}={\tt feedback}(K,\Phi)$, which when inverted gives $K={\tt feedback}(\widetilde{K},-\Phi)$. The overall controller is then $K^{*}=K_{0}+K$.

Construction of $\widetilde{K}$ uses systune where we use pole placement via TuningGoal.Poles imposing that closed-loop poles have a minimum decay of 0.9, minimum damping of 0.9, and a maximum frequency of 4.0. The controller structure is chosen as static, so that $\mathbf{x}\in\mathbb{R}^{3}$. A simulation with $K^{*}$ is shown in Fig. 9 (right) and some acceleration over the simulation for backstepping controller (left) from the same nonzero initial value is observed.

We recall that the controllers obtained in this section stabilize the system exponentially, as follows from the scheme on the right of Fig. 8, where $\widetilde{G}$, $\widetilde{K}$ are finite-dimensional and $\Phi$ is stable. This shows that $G$ is exponentially stabilizable and detectable, so that every $H_{\infty}$-stabilizing controller is also exponentially stabilizing. In particular, this applies retro-actively also to the controllers (32) and (33).

In more detail, we have to relate the trajectories of the system on the right of Fig. 8 to the trajectories on the left. Here we can follow [23], using their formulas (13), (14), to the extent that exponential decay of trajectories on the right of Fig. 8 as assured by the finite-dimensional stabilization achieved via systune leads to exponential decay of trajectories on the left of Fig. 8, which uses of course exponential stability of the infinite-dimensional part $\Phi$. While in [23] work with strictly proper $G,K$, it suffices for our present argument to suppose that all loops are well-posed. This is for instance guaranteed for a proper $K$, $\widetilde{K}$, since $\Phi$ is strictly proper.

7.3. Performance with finite-dimensional control

In this section, we show that the anti-stable wave equation (7) may be regulated satisfactorily with a simple 3rd-order finite-dimensional controller. We initialize our procedure with the controller $K_{0}$ in (32) obtained via the quasipolynomial test. Then we write the desired structure $K(\mathbf{x})$ as $K(\mathbf{x})=K_{0}+K_{1}(\mathbf{x})$, where $K_{1}(\mathbf{x})=\left[{n_{1}}/{d}\;{n_{2}}/{d}\;{n_{3}}/{d}\right]$ and deg$(n_{i})\leq 2$, deg$(d)=2$, which requires $11$ variables. This is a subclass of the class of 3rd-order controllers.

According to section 3.1, we consider the pre-stabilized system $G_{0}=G(I+K_{0}G)^{-1}$, build the closed loop ${\tt feedback}(G_{0},K_{1}(\mathbf{x}))$, and find an initial $\mathbf{x}_{0}\in\mathbb{R}^{11}$ such that $K_{1}(\mathbf{x}_{0})$ is stable and $\|K_{1}(\mathbf{x}_{0})\|_{\infty}<1/\|G_{0}\|_{\infty}$, so that by the Small Gain Theorem the loop $T(G_{0},K_{1}(\mathbf{x}_{0}))$ is stable. This is achieved e.g. by $K_{1}(\mathbf{x}_{0})=n_{0}/d_{0}[1\;1\;1]$ with $n_{0}(s)=0.3218s+0.0643$, $d_{0}(s)=s^{2}+100.1s+10$. Since $G_{0}$ has zero unstable poles, and since $K_{1}(\mathbf{x})$ is not allowed unstable poles, the Nyquist curve $1+G_{0}K_{1}(\mathbf{x})$ turns now zero times around the origin, and this is maintained during optimization.

We now use the mixed sensitivity approach of section 6.2 again, but under the form

[TABLE]

where we still have to choose the filters. The $3\times 3$ filter $W_{e}$ is chosen diagonal

[TABLE]

where the first entry is a typical low-pass, which corresponds to the output $e_{1}$. The transfer $G_{2}$ is non-minimum phase with unstable zeros at the positions $-\log(1/2)/2+jk\pi$, $k\in\mathbb{Z}$, which makes the choice of the second filter diagonal element challenging. The above choice turns out to be a good solution, as it forces tight control in the high-frequency range beyond the first unstable zero at $-\log(1/2)/2=0.346$. As third weight we choose a simple static gain $0.01$. The static filter $W_{u}=0.01$ serves to avoid unrealistic control signals.

The final controller based on (37) obtained is

[TABLE]

The final $H_{\infty}$-norm in (37) was $1.99$, with approximately $1000$ frequencies for both stability and performance. In (37), we have also constrained the controller to have a minimum decay rate of $1e-3$ and minimum damping of $0.1$ to keep control on the frequency inter-sample behavior [12]. Furthermore the constraint $|(1+K_{1}(\mathbf{x})G_{0})^{-1}|\leq 1/0.5$ stands for a disk margin of $0.5$ hence prohibiting any change in the winding number.

Simulations are shown in Fig. 10. Top left shows simulation with $K_{0}+K_{1}(\mathbf{x}_{0})$, bottom left shows the optimized controller $K_{0}+K_{1}(\mathbf{x}^{\sharp})$, achieving faster convergence and a much smaller smaller steady-state error beyond $4$ sec. Simulations of the slices $\xi=0$, $\xi=1$ and the control signal are displayed in Fig. 11 from top to bottom and confirm the previous analysis.

7.4. Gain-scheduling control

Our last study concerns the case where the parameter $q\geq 0$ is uncertain or allowed to vary in time with sufficiently slow variations as discussed in [24]. We assume that a nominal value $q_{0}>0$ and an uncertain interval $[\underline{q},\overline{q}]$ with $q_{0}\in(\underline{q},\overline{q})$ are given. The authors of [3] schedule their controller (33) using an adaptive control scheme, where the scheduling function uses a nonlinear dynamic estimate $\widehat{q}(t)\in[\underline{q},\overline{q}]$ of the anti-damping parameter.

Based on the approach in section 7.2 the following scheduling scenarios are possible. (a) Computing a nominal controller $\widetilde{K}$ at $q_{0}$ as before, and scheduling through $\Phi(q)$, which depends explicitly on $q$, so that $K^{(1)}(q)=K_{0}+{\tt feedback}(\widetilde{K},-\Phi(q))$. (b) Computing a $\widetilde{K}(q)$ which depends already on $q$, and using $K^{(2)}(q)=K_{0}+{\tt feedback}(\widetilde{K}(q),-\Phi(q))$. (c) Computing a robust controller $K_{\rm rob}$ for the entire interval.

While (a) is directly based on (6) in its finite-dimensional version based on [14, 18], see also [25], as available in systune, leading to $K^{(1)}(q)$, we show that one can also apply our approach to case (b). We use the reduction of section 7.2, see Fig. 8, to work in the finite-dimensional system $(\widetilde{G}(q),\widetilde{K}(q))$, where we now have in addition dependency on $q$, addressed by a parameter-varying design.

For that we have to decide on a parametric form of the controller $\widetilde{K}(q)$, which we chose here as

[TABLE]

and where we adopted the simple static form $\widetilde{K}_{1}(\mathbf{x})=[\mathbf{x}_{1}\;\mathbf{x}_{2}\;\mathbf{x}_{3}],\widetilde{K}_{2}=[\mathbf{x}_{4}\;\mathbf{x}_{5}\;\mathbf{x}_{6}]$, featuring a total of 6 tunable parameters. The nominal $\widetilde{K}(q_{0})$ is obtained via the synthesis technique in section 7.2. For $q_{0}=3$ this leads to $\widetilde{K}(q_{0})=[-1.049\;-1.049\;-0.05402]$, obtained via systune as in section 7.2.

With the parametric form $\widetilde{K}(q,\mathbf{x})$ fixed, we now use again the feedback system $(\widetilde{G}(q),\widetilde{K}(q))$ in Fig. 8 and design a parametric robust controller using the method of [26], which is implemented in the systune package and used by default if an uncertain closed-loop is entered. The tuning goals are chosen as constraints on closed-loo poles including minimum decay of 0.7, minimum damping of 0.9, with maximum frequency $2$. The controller obtained is (with $q_{0}=3$)

[TABLE]

with numerical values $\widetilde{K}_{1}=[-0.1102,-0.1102,-0.1053]$, $\widetilde{K}_{2}=[0.03901,0.03901,0.02855]$, and we retrieve the final parameter varying controller for the system $G(q)$ as

[TABLE]

The methods are compared in simulation in Figs. 12, 13, 14. Comparison of the simulations in Figs. 12, 13, and 14 indicates that the last controller $K_{3}(q)$ achieves the best performance for frozen-in-time values $q\in[2,4]$.

In conclusion, the study of the hyperbolic system (7) shows that optimization based on the infinite-dimensional program (6) is required to synthesize finite-dimensional controllers for (7), while its finite-dimensional counterpart based on [14] and implemented in systune is sufficient to synthesize infinite-dimensional controllers of the structure covered by Fig. 8. The major difference with parabolic systems or first-order hyperbolic systems (see e.g. [10]) is that preliminary structured stabilization, based on a suitable heuristic, cannot be verified using the Nyquist test. A very first stabilizing controller has to be found by way of some other means, but once this is achieved, the Nyquist test can be brought back to serve to control stability of the loop during optimization.