A Note on H_{2} Optimal Control Problems with Passivity Constraints ^{∗}

John E. McCarthy^{†} Clas A. Jacobson^{‡}

January 9, 1996

Abstract

This note considers the problem of determining the closest function that has positive real part on the unit disk to a given function in H_{2}. An explicit characterization of the projection is given. Remarks on potential applications using this projection result to control synthesis and system identification problems are made.

1 Introduction

The purpose of this note is to present a result connected with passivity that arises in several optimization problems associated with linear shift invariant systems. The notion of passivity has been found to be of considerable utility in capturing features of physical problems where there are natural notions of energy storage and dissipation. Of special interest are those kinds of energy notions that have been abstracted and applied to the study of the stability of feedback systems. For background in passive systems in various areas a representative sample is [Gui57, DK69] for circuit theory, [PS93] for stochastic processes, [Wil72a, Wil72b, TW91, DV75] for general passivity applications in systems theory and [Kha92, ABJ^{+}86] for specific applications of passivity in adaptive control and identification.

∗JEMcC was partially supported by the National Science Foundation grant DMS 9301508. CAJ was partially supported by the National Science Foundation grant ECS 9303356.

†Department of Mathematics, Washington University, St. Louis, Missouri 63130.

mccarthy@math.wustl.edu.

‡United Technologies Research Center, 411 Silver Lane, East Hartford, CT 06018.

jacobsca@utrc.utc.com.

The main technical result contained in this note concerns the best approximation of stable elements by passive elements. This result clarifies portions of synthesis and system identification problems that could, in principle, be solved by applications of convex programming [BB91] by specifically illustrating the structure of the projection of H_{2} onto the positive cone formed by passive systems.

We point out that, after this paper was submitted for publication, the authors learned that our main result, Theorems (2.6) and (2.7), was obtained earlier by S.D. Fisher and C.A. Micchelli in an unpublished technical report [FM86]. Notation. Let σ denote normalized Lebesgue measure on the unit circle 𝕋 in the complex plane ℂ. The Hardy space H_{2} is the subspace of L_{2}(σ) consisting of functions whose negative Fourier coefficients vanish. A function f in H_{2} can therefore be represented as a Fourier series

f(e^{it}) ∼ ∑_{n=0}^{∞} c_{n}e^{int}

where the right-hand side converges in L_{2}(σ) (and, indeed, almost everywhere on the unit circle). The function f then has a natural extension to a holomorphic function on the unit disk 𝔻, viz. f(z) = ∑_{n=0}^{∞} c_{n}z^{n}. (Following common practice the function on the circle and its extension to the disk are not distinguished typographically). For E a subset of the circle, χ_{E} denotes the indicator function of E that is 1 on E and 0 off E. For any real valued function F, the function F_{+} is max(F, 0).

Let P denote the cone in H_{2} of functions whose real parts are positive on the disk, i.e.

P = {f ∈ H_{2} : ℜf(z) ≥ 0 ∀ z ∈ 𝔻}.

Let H_{2,ρ} denote the Hardy space H_{2} on the disk of radius ρ and P_{ρ} the associated cone of positive real functions. Let H_{2,ρ}(M) denote the ball of elements satisfying ∥f∥_{2,ρ} ≤ M. Let H_{∞}^{R} denote the subset of H_{∞} that has real values for real arguments (conjugate symmetric functions).

2 Main Result

The purpose of this section is to find, for a given function f in H_{2}, the best approximant in P. Note that because P is a closed and convex subset of a Hilbert space every function will have a best approximant in P and this approximant will be unique [Lue69]. We shall use the notation that f is the function to be approximated and u is the best approximant in P.

Note that ℑ(f(0)) = ℑ(u(0)), so without loss of generality ℑ(f(0)) = 0. Every function in H_{2} can be written as the sum of its real and imaginary parts, and the norms of these are related by

∫ (ℜ[f − f(0)])^{2}dσ = ∫ (ℑ[f − f(0)])^{2}dσ.

So

∥f − u∥^{2} = ∫ (ℜ[f − u])^{2} + (ℑ[f − u])^{2}dσ
= ∫ (ℜ[f − u])^{2} + (ℜ[f − u − (f(0) − u(0))])^{2}dσ
= 2 ∫ (ℜ[f − u])^{2}dσ − (f(0) − u(0))^{2}

Let F denote ℜ(f) and U denote ℜ(u). The restriction that ℜ(u) be positive on 𝔻 is equivalent to requiring U to be positive on the circle. Finally note that

f(0) − u(0) = ∫ f − u dσ = ∫ F − U dσ.
□

Combining the above the projection problem is then equivalent to the following one:

Given a real-valued function F in L_{2}(σ), find the positive function U that solves the extremal problem of minimizing

2 ∫ (F − V)^{2}dσ − (∫ F − V dσ)^{2}  (2.1)

over all positive functions V.

Lemma 2.2 With notation as above, U − F ≥ 0 almost everywhere.

Proof:Let E = {t ∈ 𝕋 : U(t) − F (t) < 0}. Let ε > 0, and let V = U + εχ_{E}. Let |E| denote σ(E). Then

2 ∫ (F − V)^{2}dσ − (∫ F − V dσ)^{2} =

2 ∫ (F − U)^{2}dσ − (∫ F − Udσ)^{2} + ε^{2}|E|

−4ε ∫_{E} F − U dσ + 2ε|E| ∫ F − Udσ
= 2 ∫ (F − U)^{2}dσ − (∫ F − Udσ)^{2} + ε^{2}|E|
−2ε(2 − |E|) ∫_{E} F − U dσ + 2ε|E| ∫_{E^{c}} F − U dσ

As U is the infimum it follows that that the sum of the first-order terms in ε be positive, so

−ε(2 − |E|) ∫_{E} F − U dσ + ε|E| ∫_{E^{c}} F − U dσ ≥ 0.  (2.3)

Therefore

|E| ∫_{E^{c}} F − U dσ ≥ (2 − |E|) ∫_{E} F − U dσ  (2.4)

But the integrand on the right-hand side of (2.4) is greater than zero, and the integrand on the left is less than or equal to zero. Therefore the measure of E must be zero, as required. □

A similar argument yields the following:

Lemma 2.5 Let Γ = {t ∈ 𝕋 : U(t) > 0}. Then

∫_{Γ} U − F dσ = \frac{|Γ|}{2} ∫_{𝕋} U − F dσ.

Proof:Let δ > 0, and let Γ_{δ} = {t ∈ 𝕋 : U(t) > δ}. The same variational argument as in the previous lemma can be applied for V = U + εχ_{Γ_{δ}}, but V is now positive for all ε > −δ. Therefore inequality (2.3) holds for both positive and negative values of ε, so the left-hand side must be zero:

|Γ_{δ}| ∫_{Γ_{δ}^{c}} F − U dσ = (2 − |Γ_{δ}|) ∫_{Γ_{δ}} F − U dσ.

Letting δ tend to zero and doing a little algebra gives the desired result. □

The main result can now be stated. For any real-valued function G, the notation G_{+} stands for max(G, 0).

Theorem 2.6 The optimal solution U is given by U = (F + c)_{+}, where c = \frac{1}{2} ∫ U − F dσ is a non-negative constant.

Proof:Consider first the case when Γ = {U > 0} has positive measure. Claim: On Γ, the function F − U is constant. Proof of claim: Deny. Then there exist ε > 0 and sets E_{1} and E_{2} in Γ with the following properties:

1. (i) |E_{1}| = |E_{2}| > 0.

2. (ii) U > ε on E_{2}.

3. (iii) inf_{t∈E_{2}}{U(t) − F (t)} > sup_{t∈E_{1}}{U(t) − F (t)} + 2ε.

Now let V = U − εχ_{E_{2}} + εχ_{E_{1}}. The function V is positive, and ∫ F − V dσ = ∫ F − U dσ. Moreover

∫ |V − F |^{2} dσ = ∫ |U − F |^{2} dσ + 2ε^{2}|E_{1}|
+2ε ∫_{E_{1}} U − F dσ − 2ε ∫_{E_{2}} U − F dσ
≤ ∫ |U − F |^{2} dσ + 2ε^{2}|E_{1}| − 4ε^{2}|E_{1}|.

This would contradict the minimality of U in (2.1), so the claim is established.

This proves that U = (F + c)_{+} for some constant c, at least on the set Γ ∪ {F ≤ −c}. If this set were not of full measure, one could choose sets E_{2} in Γ and E_{1} in {−c < F ≤ 0} ∩ {U = 0} that satisfy properties (i) through (iii), and the same computation as above would yield a contradiction.

Therefore U = (F + c)_{+}. By Lemma 2.2 the constant c is non-negative, and by Lemma 2.5 it is \frac{1}{2} ∫_{𝕋} U − F dσ.

Now consider the case where |Γ| = 0. This can only happen when F is non-positive (though this is not sufficient to make U identically zero). To finish the proof of the theorem, we have to show that in this case \frac{1}{2} ∫ F dσ is greater than or equal to the essential supremum of F - call this number L (when |Γ| = 0, the choice of c is not unique, in general; we are claiming that the choice −\frac{1}{2} ∫ F dσ will work).

For δ > 0, let E = {|F − L| < δ}, and let V = εχ_{E}. The minimality of U = 0 in (2.1) gives, for ε positive and tending to zero, that

2|E| ∫ F ≥ 4 ∫_{E} F dσ.

Letting δ tend to zero gives that ∫ F dσ ≥ 2L, as desired. □

Theorem 2.6 only gives c in terms of the (still unknown) function U. The next theorem addresses how to find c directly; note that it works whether or not |Γ| = 0, and that once one finds c one knows U from the previous theorem.

Theorem 2.7 The number c = \frac{1}{2} ∫ U − F dσ is the unique point at which the decreasing function

G(x) = x|{F > −x}| − ∫_{F ≤−x} F dσ − 2x,

defined for x ≥ 0, changes sign.

Proof:To show that G is decreasing, consider, for positive h,

G(x + h) − G(x) = (x + h)|{F > −x − h} − ∫_{{F ≤−x−h}} F dσ − 2x − 2h
−x|{F > −x}| + ∫_{{F ≤−x}} F dσ + 2x
= −x|{−x − h < F ≤ x}| + ∫_{{−x−h<F ≤x}} F dσ
+h|{F > −x − h}| − 2h
≤ −2x|{−x − h < F ≤ x}| − h
< 0,

so G is decreasing (the above also shows that G is lower semi-continuous, with jumps at precisely those points x for which the level set {F = −x} has positive measure). G(0) is greater than zero, unless the original function F is non-negative and therefore equal to U. As x tends to ∞, G(x) + x is o(x), so G(x) tends to −∞. Therefore there is a unique point at which G(x) changes sign.

To see that this point is c, consider Lemma 2.5, which can be rewritten as

\frac{1}{2 − |Γ|} ∫_{Γ^{c}} U − F dσ = \frac{1}{2} ∫_{𝕋} U − F dσ.

The right-hand side of this equation equals c, and, as Γ = {U > 0} = {F > −c}, it follows that

\frac{1}{2 − |\{F > −c\}|} ∫_{\{F ≤−c\}} −F dσ = c.

Cross-multiplying gives G(c) = 0. □

Example. Consider the function f(z) = z. Parameterize the unit circle by the interval [0, 1], so F (t) = cos(2πt). By Theorem 2.7, c is the solution of the equation

c (1 − \frac{1}{π} cos^{−1}(−c)) + \frac{1}{π} \sqrt{1 − c^{2}} − 2c = 0

This can be solved numerically, and one gets that c = 0.217234 … Therefore the positive function U that minimizes 2.1 is U(t) = (cos(t) + .217234 …)_{+}.

To solve the original problem, one must compute the harmonic conjugate Ũ of U, and then u = U + iŨ. One way to do this is in terms of Fourier series: if the function u(z) = ∑_{n=0}^{∞} c_{n}z^{n}, then c_{0} = ∫_{0}^{1} U(t)dt, and for n ≥ 1,

c_{n} = 2 ∫_{0}^{1} U(t)e^{−2πint}dt

In this particular case, the first few terms in the Fourier series are (to five decimal places):

u(z) = .43447 + .6372z + .19736z^{2} − .04287z^{3} − .02830z^{4} + .02249z^{5}

+.00515z^{6} − .01292z^{7} + .00151z^{8} + .00723z^{9} − .00353z^{10} + …

Note that if one approximates u not by the partial sums of its Fourier series, but by their Cesàro means, then the approximants themselves will all have positive real part, by the positivity of the Fejér kernel. The function F (t) = cos(2πt), the best approximant U, the 8^{th} order partial sum of the Fourier series S and the 8^{th} order Cesàro mean C are shown in Figure 1.

Figure 1: Best Approximation using Positive Real Elements

Simple Functions. Assume that F is a finite linear combination of characteristic functions,

F = ∑_{i=1}^{N} α_{i}χ_{E_{i}}  (2.8)

where the E_{i}’s are pairwise disjoint measurable sets whose union is the whole circle, and α_{1} > α_{2} > … > α_{N}. Let t_{i} = |E_{i}|. For r ≤ N, consider the system of linear equations

( \begin{matrix}1 − \frac{1}{2}t_{1} & −\frac{1}{2}t_{2} & … & −\frac{1}{2}t_{r} \\ −\frac{1}{2}t_{1} & 1 − \frac{1}{2}t_{2} & … & −\frac{1}{2}t_{r} \\ ⋮ & ⋮ & ... & ⋮ \\ −\frac{1}{2}t_{1} & −\frac{1}{2}t_{2} & … & 1 − \frac{1}{2}t_{r} \end{matrix} ) ( a_{1} ; a_{2} ; ⋮ ; a_{r} ) = ( α_{1} − \frac{1}{2} ∑_{j=1}^{N} t_{j}α_{j} ; α_{2} − \frac{1}{2} ∑_{j=1}^{N} t_{j}α_{j} ; ⋮ ; α_{r} − \frac{1}{2} ∑_{j=1}^{N} t_{j}α_{j} )

Let r be the largest integer for which these equations have a solution a_{1}, …, a_{r} for which all the a_{i}’s are positive. Then the best approximant to F is

U = ∑_{i=1}^{r} a_{i}χ_{E_{i}}.

Proof:The best approximant must be of the form

U = ∑_{i=1}^{N} b_{i}^{2}χ_{E_{i}}

where b_{i}^{2} = (α_{i} + c)_{+}, for some c. Taking partial derivatives of

2 ∫ (F − U)^{2}dσ − (∫ F − Udσ)^{2}  (2.9)

with respect to b_{i} and setting equal to zero gives either b_{i} = 0 or

(b_{i}^{2} − α_{i})t_{i} − \frac{1}{2}t_{i} [∑_{j=1}^{N}(b_{j}^{2} − α_{j})t_{j}] = 0  (2.10)

The restriction that b_{i}^{2} ≥ 0 means that (2.10) cannot, in general, be satisfied by all the b_{i}’s. Now suppose r is the largest index for which (2.10) can be solved for all i ≤ r, setting b_{i} = 0 for i > r. Let a_{i} = b_{i}^{2}. Then ∑ a_{i}χ_{E_{i}} must give the optimal solution. For there cannot be more non-zero terms in the optimal solution, by construction. And there cannot be less either, as the a_{i}’s give the unique critical point of (2.9) among all linear combinations of χ_{E_{1}}, …, χ_{E_{r}} without the restriction that the coefficients be positive. □

This method, since it only involves solving linear systems, is an efficient way of finding U if F is a simple function (i.e. a function of the form (2.8)).

Note that one candidate for u is v = F_{+} + iF̃_{+}. The results above show that this is not in general optimal, but it is not too bad: it can be shown that ∥f − v∥ ≤ \sqrt{2}∥f − u∥.

It may also be useful to solve a related problem, namely to minimize (2.1) over all functions V that are positive not necessarily on the whole circle, but just on a subset^{1} E. The solution can be found using similar ideas to those utilized above; we shall not write out the details of the proof, but just give the result here.

Theorem 2.11 Let F be a real-valued function in L_{2}(σ), let E be a set of positive measure in 𝕋, and let U be the function that minimizes

2 ∫ (F − V)^{2}dσ − (∫ F − V dσ)^{2}

over all functions V that are positive on E. Let

λ = \frac{2 + |E^{c}|}{2 − |E^{c}|}.

1This problem arises for example in the situation where a passive element is cascaded with a delay, i.e. in continuous time this would happen with a transfer function of the form P (s)e^{−sT}, that is common in (practical) collocated systems. This projection result is relevant for system identification of such systems as considered in section 3.2.

Then

U = \{ \begin{matrix}(F + c_{1})_{+} & on E \\ F + c_{2} & on E^{c} \end{matrix}

where c_{1} is the unique zero of the function

G(x) := λx|{F > −x}| − λ ∫_{F ≤−x} F dσ − 2x,

and

c_{2} = \frac{∫_{E} U − F dσ}{2 − |E^{c}|}
= \frac{c_{1}|\{F > −c_{1}\} ∩ E| + ∫_{\{F ≤−c_{1}\}∩E}(−F)dσ}{2 − |E^{c}|}

3 Applications

3.1 Stabilization by Stable Compensators

The purpose of this section is to point out how the approximation result given in section 2 can be utilized to parametrize a set of stable stabilizing compensators for a quite general (and physically motivated) control configuration.

The discussion here is on a simple control configuration (which arises in applications, see for example [NE92]). Generalized results based on this idea are certainly possible and are left as future investigations. Consider the control configuration shown in Figure 2.

For simplicity assume that all elements (P, C, F) are stable (that is, in H_{2}) and F is passive, that is, F ∈ P, then the control problem is to find W such that the system is stable and the map r ↦ y is minimized (that is, an H_{2} model matching problem of the form considered in [Fra82]).

Either directly or by converting this to the standard problem configuration [Fra87] the set of all stabilizing W can be given by the so called Q-parameterization as

Q = \frac{W}{1 − W F}

and the map r ↦ y is given by

QC − P

[figure: F; W; C; P; +; +; +; +]

Figure 2: Control Configuration for Model Matching Problem

Consequently the performance as measured by the H_{2} norms of r ↦ y, ∥QC − P ∥_{2}, may be optimized by an inner-outer factorization of C followed by a standard L_{2}/H_{2} approximation problem; moreover, since P ⊆ H_{2} the approximation problem fits into the framework of section 2 if the parameter Q is desired to be chosen to be passive.

The point of this discussion is exactly the consequences of choosing Q to be passive: if Q ∈ P then the controller W will be passive and hence stable. This follows immediately from the control configuration under consideration and the standard passivity theorem [DV75, p. 182]. This observation certainly has implications for robustness (when in particular F may be perturbed by passive perturbations) which are left for future investigations; however, what should be especially noted is the simplicity of this type of stable stabilizing controller result as contrasted with the general case addressed in [GP89] and the algebraic case (that is, no optimal performance) of just stable stabilization in [YBL74]. This can be summarized as follows.

Theorem 3.1 Consider the control configuration shown in Figure 2. Assume that P, C ∈ H_{2}(𝔻), F ∈ P. The H_{2} optimal solution of the model matching problem

min_{Q∈P} ∥CQ − P ∥_{2}

will result in a stable passive compensator

W = \frac{Q}{1 + F Q}

It should be noted that the simplifications of this result as contrasted with others is due to the strength of the following assumptions: (1) the control configuration and (2) the passivity of F. Note also that the results of the previous section apply directly only if C is inner - this corresponds to f being the projection onto H_{2} of C̄P. If C has a non-constant outer factor, then one must solve the problem of finding the projection onto the passive elements in a weighted H_{2} space.

3.2 System Identification of Passive Systems

The complete problem setting that we consider in this section is contained in [HJN91, HJN93] (see [TWW88, GW59] for a detailed discussion of the general situation of optimal recovery). A summary of the problem statement is as follows. Given a function h ∈ H_{∞}(𝔻_{ρ}) with ρ > 1, ∥h∥_{∞,ρ} ≤ M and a set of corrupted frequency samples E_{n}(h, ϵ) = {h(e^{iω_{k}}) + η_{k}, k = 1, …, n, η ∈ ℂ^{n}, ∥η∥_{∞} ≤ ϵ} and ω_{k} selected as the n roots of unity. Construct an approximation by an operator (in general nonlinear) A : ℝ_{+}×ℝ_{+}×ℂ^{n}×ℝ_{+} → H_{∞}(𝔻) by

f = A(M, ρ, ϵ, E_{n})

and an a priori error bound

∥h − f∥_{∞}

suitable for use in a subsequent design procedure.

Space does not permit us to discuss the solution at length. However we wish to point out that in solving the H_{2} version of the problem, it is possible to take the positive real character of a function to be identified — the cone property — into algorithms as easily as the usual subspace constraints. Indeed, one can first find the H_{2} solution, and then project it onto the cone of passive elements (or the cone of functions whose real parts are positive on some fixed arc). This becomes more difficult when working in the H_{∞} case, where the projection is not well understood (and indeed has no unique solution in general). But, and this is the reason for introducing 𝔻_{ρ}, the H_{∞}(𝔻) norm is dominated by the H_{2}(𝔻_{ρ}) norm, so this can be used to obtain a good (if not necessarily optimal) solution. See [MU88] and [HJN93].

Acknowledgment

The authors are grateful to Alexander Pruss for some useful comments.

References

[ABJ^{+}86] B.D.O. Anderson, R.R. Bitmead, C.R. Johnson, P.V. Kokotovic, R.L. Kosut, I.M.Y. Mareels, L. Praly, and B.D. Riedle. Stability of Adaptive Systems: Passivity and Averaging Analysis. MIT Press, Cambridge, MA, 1986.

[BB91] S.P. Boyd and C.H. Barratt. Linear Controller Design. Prentice-Hall, New York, 1991.

[DK69] C.A. Desoer and E.S. Kuh. Basic Circuit Theory. McGraw-Hill, New York, 1969.

[DV75] C. A. Desoer and M. Vidyasagar. Feedback Systems: Input–Output Properties. Academic Press, New York, 1975.

[FM86] S.D. Fisher and C.A. Micchelli Minimal norm interpolation in H^{2} with nonnegative real part. IBM Technical Report RC 12351 (# 55390), December 1986.

[Fra82] B.A. Francis. On the Wiener-Hopf approach to optimal feedback design. Systems and Control Letters, 2:197–201, December 1982.

[Fra87] B. A. Francis. A Course in H_{∞} Control Theory. Springer–Verlag, 1987.

[GP89] C. Ganesh and J.B. Pearson. H^{2}-optimization with stable controllers. Automatica, 25(4):629–634, 1989.

[Gui57] E.A. Guillemin. Synthesis of Passive Networks. Wiley, New York, 1957.

[GW59] M. Golomb and H.F. Weinberger. Optimal approximation and error bounds. In R.E. Langer, editor, On Numerical Approximation, pages 117–190. Univ. of Wisconsin Press, Madison, 1959.

[HJN91] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. Control-oriented system identification: A worst-case/deterministic approach in H_{∞}. IEEE Transactions on Automatic Control, 36(10):1163–1176, October 1991.

[HJN93] A.J. Helmicki, C.A. Jacobson, and C.N. Nett. Identification in H_{∞}: Linear algorithms. IEEE Transactions on Automatic Control, 1993.

[Kha92] H. K. Khalil. Nonlinear Systems. Macmillan, New York, 1992.

[Lue69] D. G. Luenberger. Optimization by Vector Space Methods. Wiley, New York, 1969.

[MU88] C.A. Micchelli and F.I. Utreras. Smoothing and interpolation in a convex subset of a Hilbert space. SIAM J. Sci. Stat. Computing, pages 728–746, 1988.

[NE92] P.A. Nelson and S.J. Elliott. Active Control of Sound. Academic Press, London, 1992.

[PS93] S.U. Pillai and T.I. Shim. Spectrum Estimation and System Identification. Springer-Verlag, New York, 1993.

[TW91] H.L. Trentelman and J.C. Willems. The dissipation inequality and the algebraic Riccati equation. In S. Bittanti, A.J. Laub, and J.C. Willems, editors, The Riccati Equation, pages 197–242. Springer-Verlag, 1991.

[TWW88] J. F. Traub, G. W. Wasilkowski, and H. Wozniakowski. Information-Based Complexity. Academic Press, 1988.

[Wil72a] J.C. Willems. Dissipative dynamical systems part I: General theory. Arch. Rational Mech. Anal., 45:321–351, 1972.

[Wil72b] J.C. Willems. Dissipative dynamical systems part II: Linear systems with quadratic supply rates. Arch. Rational Mech. Anal., 45:352–393, 1972.

[YBL74] D.C. Youla, J.J. Bongiorno, and C.N. Lu. Single-loop feedbackstabilization of linear multivariable dyanmical plants. Automatica, 10, 1974.