Data-based synthesis of a multivariable linear-quadratic regulator

Auromoricn, Vol. 32, No. 3, pp. 403407,

Pergamon

1%

Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved OOOS10!28/% $15.00 + 0.00

00051098(95)00142-5

Brief Paper

Data-based Synthesis of a Multivariable LinearQuadratic Regulator* JENQ-TZONG Key Words-Linear-quadratic

regulator; computer-aided

1. Introduction The linear-quadratic (LQ) optimal regulator is commonly designed by means of a state feedback law, derived as follows. Consider a linear sysystem with an fi X 1 output vector y, an m X 1 input vector II and an n X 1 state vector x: y(k) = Cx(k),

design; time-domain data analysis.

duplicate the function of a regulator designed by the state feedback technique. As a result, the need for state estimation is eliminated. In the proposed method, LQ design is based purely on input and output data. The control input sequence and the output sequence of an optimal system are computed numerically without explicit identification of a parameter&d plant model; a correlation equation between the LQ response data and the proposed output feedback controller makes it possible to compute the controller parameters when the plant model has not yet been identified. As a result, plant model estimation and LQ controller synthesis are combined into a single design procedure.

Abstract-A new method for the design of a linearquadratic optimal regulator directly from input-output data is proposed. The method can be used to design a multivariable regulator when the plant is open-loop stable, with more outputs than inputs, and is left-invertible. The major advantages of the method are threefold. First, it involves an output feedback control law; hence, no state estimation is needed for implementation. Second, the computation of this controller is possible without explicit identification of the plant model parameter. Third, a reduced-order controller, if one exists, may be synthesized without additional effort to reduce the plant model.

x(k + 1) = Ax(k) + h(k),

H. CHANt

2. Optimal output feedback regulator Denoting the z transforms of {x(k)}, {u(k)} and {r(k)) by X(z), ii(z) and r(z) respectively, we may write the state feedback control law for the z domain as

(1)

ii(z) = -fi(z)

where A, B and C are constant matrices of appropriate dimensions. The LQ design calls for minimization, with respect to u, of the cost functional

From (1) we infer H,,(z)ii(z), where

that

Z(z) =H,(z)ii(z)

H,(z) = (al -A)-‘& J = t $ [y(/~)~Qy(k) + uTRu(k)], k-0

Q 20,

R >O.

+?(z).

(3) and

f(z)=

H,,(z) = C(z1 - A)-‘B.

(2) In this paper, we assume that a left inverse of H&z) exists. If we denote this left inverse by HiL(z), the output feedback equivalent of (3) may be written as

As p-+ m, a steady-state feedback control law u(k) = -Kx(k) + r(k) results, where K is a 1 X n gain vector and r is an m X 1 command vector. This method, however, has several shortcomings. First, implementation of the state feedback control law requires a state estimator (Sage and White, 1977), the synthesis of which requires additional effort. Also, the inclusion of a state estimator in the system diminishes the robustness of the controller (Doyle and Stein, 1979). Second, the computation of the control law requires a fairly precise knowledge of the plant model. Often, this knowledge is obtained from plant input and output data through system identification (Ljung and Siiderstriim, 1983; Astrom and Wittenmark, 1989). System identification, however, adds complexity to the task of control system synthesis. Multivariable system control using output feedback to achieve various closed-loop system designs is possible. In this paper, we show that, for left-invertible systems, an LQ regulator designed using an output feedback control can

ii(z) = -F(z)y(z)

+ T(z),

F(Z) = KH,(z)GL(z),

(4)

where F(z) represents an output feedback controller. We then express H,(z) in coprime fraction form as H,(z) = %&)/a(z), where a(z) is a polynomial and 9(z) is an n x m polynomial matrix. Then HyL(z) can be expressed as a(~)%;~(z)/q(z), where q(z) is a polynomial and %iL(z) is an m x fi polynomial matrix. As a result, F(z) can be expressed as F(z)

=

z,

P(z) = KYe,(z)%GL(z).

For F(z) to be implementable, the rational function P(z)/g(z) has to be proper. Also, the closed-loop equation of the output feedback system can be derived from (4) and (5) as follows: Y(z) = ]a(z)rr(z)l+

*Received 29 November 1993; revised 30 August 1994, revised 11 July 1995; received in final form 22 August 1995. This namer was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor A. Tits under the direction of Editor T. Basar. Corresponding author Professor Jenq-Tzong Chan. E-mail 7575 ext. 63681; Tel. +886 6 275 pub@alex,iaa,ncku,edu,tw. t Institute of Aeronautics and Astronautics, National Cheng Kung University, Tainan, Taiwan, 70101, Republic of China.

~,l~z~~~z~l-‘sB,,~z~B~z~~~z~. (6)

On the other hand, the closed-loop state feedback design is j(z) = i&(z)f(z),

system in the optimal

H,(z) = C(zl -A + BK)-‘B.

(7)

By comparing (6) and (7), we obtain the following equation: ]o(z)cr(z)l+

%(z)~(z)1-‘%(z)g(z)

= H,(z).

(8)

From (8) we see that a transmission zero, g(z), which does 403

Brief Papers

404

not exist in the plant equation, may appear in the closed-loop dynamics of a state feedback system, a situation that has been proven impossible in linear systems (Chen, 1984). Hence, cancellation of g(z) must occur in (6). If the system has an equal number of inputs and outputs then &TL(z) = &‘(z), and q(z) will be unique and specific. In this case, unstable pole-zero cancellation will result, and the output feedback controller cannot be used, if g(z) is non-Hunvitz. It has been shown (Chan. 1995) that if (i) ti = m + 1, (ii) CB is left-invertible and (iii) H,,(z) has no transmission zeros (Chen, 1984) than an HdL(z) with an arbitrary polynomial in the denominator always exists. In such cases, g(z) may be freely assigned, and the problem of unstable pole-zero cancellaton is eliminated. Moreover. a proper S(z)/g(z) will be ensured if the matrix product CB is left-invertible. If these conditions are satisfied, a workable output feedback controller always exists. We may state this output feedback controller as follows:

ii(z) =

-q(z) dz)

+r(z).

Y(z) =

if;zI /=(I

‘,

(9)

where e are m x Si constant matrices. In the following, we present a numerical process for the computation of F; from input-output data without the explicit identification of H,,(z ). 3. Details of the design problem

formulation for finite-time LQ 3.1. Numerical minimization. Equation (2) can be expressed in matrix form

in terms of input and output data:

(10) where Q and R are block-diagonal matrices with diagonal blocks equal to Q and R respectively. Note that the elements of U start at u(O) and end at u(p - 1). while those of Y start at y(1) and end at y(p). This is due to causality within the plant dynamics, which dictates that u(k) cannot affect y(l) for all 15 k (Franklin and Powell, 1980). In general, the output of a controlled system consists of the input response plus the initial state response caused by the non-relaxed initial condition. Therefore, if we denote the unit pulse response sequence of I&,(z) by {H(k)}t and the initial state response sequence as {a(k)} then the following equation is inferred:

As a result, the performance index (10) becomes J = :(WU + h)‘Q(WU

+ 72) + :U’IwUJ.

(12)

The optimal regulating input is derived by minimizing J with respect to U: U = -(W’l‘QW + W) ‘W’K@.

(13)

It is also possible to replace H with some multicolumn matrix,

Accordingly. matrix.

111will also be replaced

by a multicolumn n,(O)

.

U- [LJ,

;

U,,], u, =

,

[ Y(P) 1 in which each m X 1 vector sequence {ui(k)} represents the optimal input sequence corresponding to Idi(k As a result, (13) may be rewritten as U = -(WrQW + lR-+H%@

(14)

Note that (14) is computable even when the plant model is not known. The computation requires only the data for {a,(k)} and {H(k)}. However, it may not be desirable to directly measure {H(k)}. The plant response to a more appropriate test command may be easier to sample. In this case, {H(k)} may be computed from the inverse convolution formula

[H(l)

“’

wP)l=[‘w)

“’

X(P)1

T(O)

T(p-1)

-’

X

1

“.

:

TiO)

WI

I

where {T(k)} is an m X m matrix sequence representing the test command and {x((k)} is an m X m matrix sequence representing the plant response. The data within {x((k)} is sampled, one column at a time, from m repeated tests, each test using a different column of {T(k)} as test input. The chosen sequence {T(k)} must ensure a well-defined matrix inversion in (15). 3.2. Steady-state opcmal solution. Consider the case in which the data within Z in (14) are generated by m repeated tests, with each test using a unit-pulse command applied to a different input channel at k = -1. In this case, a)(&) will be the ith column of H(k + 1) for all k SO and Z in (14) becomes

Then. denoting U(O)

;

uLi=

) where U(1) are m X m matrices,

[ WP - 1) I the data sequence {U(k)} computed

from (14) becomes the optimal control input response to a unit pulse command. The LQ formulation guarantees the stability of the optimal design (Sage and White, 1977). Hence the convergence to zero of {U(k)} is ensured. A steady-state solution is also guaranteed if the final time of the optimization is sufficiently large (Franklin and Powell, 1980). For sufficiently large p, it follows that for any small real number p > 0, there always exist finite positive integers p and k such that /(U(k) 11 < p for all p 2 k > k.$ In this case, the data sequence {U(k)} is a steady-state optimal control input response to a unit pulse command. Two approaches may be used for the computation of the steady-state solution. (i) One can follow the procedure of Franklin and Powell (1980) for model-based LQ design. I. Select

an appropriate

accuracy

threshold

p and a proper

value for p. 2. From W from (11). 2 from (16). and compute fi from (14).

where the r?i X 1 vector sequences {d,(k)} are m sets of initial state respones, corresponding to different initial conditions.

3. Check to see if IiC/(k) 1)< p for all k 2 l$p. If this condition is not met, double the value of p and go to step 2. (In this procedure, the value of k is made equal to lg. However, it may be made equal to any other fraction of p as well.) (ii) When, at a sufficiently small value of p, an appropriate

t {H(k)} is an ri7 X m matrix sequence. unit-pulse responses of the plant.

*The norm the magnitudes

Its elements

are the

of a matrix M, l/MiI, is defined of its elements.

as the sum of

Brief Papers value of k has been determined (see Remark 2 below), we may regard U(k) = 0 for all k > k. Then, using a formula derived in the Appendix, a recursive estimation of {U(k)} may be performed as follows:

w (j= ;WV ,1

sl,+, = aj, - F,+I@;+,~~+,

405

and, because g(z) = zf_.,, gizemi, D = [D(O) D(k) = c g@(k i=O

(17)

8(k) = {;

where 4 is an estimate, at time p, of the data vector

Q[WP + 2) + @p+&A H(p -k)].

Remark 1. Since (17) is a recursive least-squares formula, its convergence is assured (Ljung, 1987). Also, the definition of FP ensures the existence of F;’ for all p ~0; hence the problem of persistent excitation (Anderson and Johnson, 1982) does not arise in (17). Remark 2. The assumption of a stable plant implies that lim,,, H(k) = 0. Hence, for sufficiently small p >O, it is always possible to determine from the sequence {H(k)} a finite integer k>O such that IIH(k)j(

c. An appropriate value for k many be a multiple of k

If {U(k)} is the optimal control input response to a unit-pulse command then the optimal unit-pulse output response of the LQ design is equal to the convoltion of the sequence {U(k)} and the unit pulse plant response. If {Y(k)} is an 3 Xm matrix sequence representing the optimal unit-pulse output response then the following equation holds: Y(k)J=[I

U(0)

x[H(l)

U(k-

1::

?],

P(z) = I

O(z),

-

ii(z) = - &Y(z) e(z)

,$.+

=y

(19)

(21)

Then the following data correlation equation is inferred from the inverse z transform of (21): PW=D,

(22)

-‘F(z).

P=[F,

...

F8],

W=

.. ‘..

(25)

(26)

’

will be invertible and will have a stable inverse. For an arbitrary value of e, let * + S(z) o(e)&(z)

= G(z) + E(z),

(27)

where E(z) denotes the error in s(z) that causes {I + [@z)/e(z)]&(z)}-’ to differ from the LQ solution. By using the generalized Nyquist stability criterion (Pate1 and Munro, 1982), we may derive the following sufficient condition for the stability of the closed-loop system: I+:

II~{E(z)G(z)-‘}Il,=~e< denotes

1,

(28)

the maximum eigenvalue of

E(z)G(z)-‘.

Furthermore, (23) by

if we denote the correlation error matrix of

. . E(k)] =

(29)

then a solution at a correct value of C will cause E = 0. For an arbi$ary value of e, let the z transform of {E(k)} be denoted by E(z). Then (20) and (21) imply that @@@[U(z)-z]+WT(z) ZP

ZP

where Y(1)

(24)

since the plant dynamics dictates that Furthermore, f(z) = II&z)ii(z), and a stable I&(z) is assumed, the stability of the design is ensured if the closed-loop dynamics of ii(z) are stable. For a correct value of e, the LQ formulation ensures that the transfer-function matrix

(20)

qz)].

’1

ii(z) = I+%&(z) [ e(z)

IE=[E(l) [I -

S(z) = 2 fiz’f’ i=o ’

with 8 being the elements of P. Then the closed-loop dynamics of the control input becomes

where i{E(z)G(z)-I}

Y(z) =

+ T(z),

G(z)=,+%&(z) (l(z)

where u(z) and Y(z) represents the z transforms of the data sequences {U(k)} and {Y(k)} respectively. An alternative expression of (20) may be written as follows:

y

(23)

3.4. The determination of L It has been shown in Section 2 that an output feedback controller that can achieve LQ design always exists for some e> 0. Since the plant model is assumed to be unknown, trial values of e, starting from e = 1, are used in (23). Consequently, if a reduced-order controller exists, it will manifest itself in the trial process. The trial process requires a stability criterion, derived as follows. Given a computed P, the proposed controller can be implemented as follow:

l)]

where the m X m identity matrix I represents the unit-pulse command applied at k = -1. 3.3. Numerical computation of 4. If {U(k)} is the control input, and {Y(k)} the output of the optimal system in response to a unit pulse command, then, according to the output feedback control law given by (9), the following correlation equation holds:

g

I: ;;;;

3. Although (23) has been derived using a least-squares technique, other more advanced techniques, such as SVD- or QR-based methods (Golub and Van Loan, 1985), may be used. Also, recursive implementation of (23) is possible. See, for instance, Ljung (1987) and Astram and Wittenmark (1989) for the detailed procedure needed to convert (23) into recursive formulas.

Note that matrix inversion is avoided in the recursive formula. Recursive estimation will also produce a consistent and bias-free result when the data within {H(k)} have been corrupted by noise (Franklin and Powell, 1980). (The noise, however, must be white and unbiased.)

...

- U(k -i)],

Remark

and

‘[Y(l)

-i)

$ = DW(@I@-‘.

[

Qp= [H(p) . . .

D(k - l)],

As a result, a least-squares estimate of P, denoted by P, may be computed as follows:

FP is a time-varying gain defined as

8 p+, =

.

=~[u(z)-I]+~Ho,~z)LI(L)

406

Brief Papers

Because the sequence {U(k)} is an optimal LQ solution to a unit-pulse command, the relation U(z) = G(z)-’ holds. By substitution into (30), we infer that g(z) = “OE(z)G(z)-I. 2’

When a noise signal is present, H,,(z) is not computable. In this case, the value of H&z) for all z on the unit circle may be estimated instead from the DFT of {H(k)}. If we thus rewrite (36) in terms of Q(z), we arrive at the equation

(31)

Consequently, (28) may be rewritten as (32) In practice, &ej@) may be approximated by the discrete Fourier transform (DFT) of the error sequence {E(k)}. Also, the Frobenius norm (Reid, 1983) of a square matrix normally presents an upper bound for its characteristic values, and is numerically much simpler to compute. If we denote the DFT of {Elk11 . . bv E(k)I for 0 < k I k. and the Frobenius norm of E(k) by l~~(k)llr. then a computable stabihty criterion for the designed system emerges: I_

_.

.

,,

(33) Note that {E(k)} and {E(k)} may be computed from (29) for each vahre of e tested. A stable solution is obtained when (33) is satisfied, at whatever value e may have. Hence it is unnecessary to guess the value of P exactly, and concerns about its uncertainty are removed. 3.5. Numerical errors in the formulation. At least two possible sources of error exist in the formulation: the noise corruption in {H(k)} and the discarding of U(k) for all k > k in the recursive computation of 111.Analysis of these errors is given as follows. First, let the error vector of (14) be defined as N=

N(I) ; I N(P) 1

+ % 3Q[Wz) + %)I + MWz) + [H,,(z) + ~(z,l’Qz~Ho,k~ + %)I,

+ fib)1 (35)

where the prime denotes the transpose of the complex conjugate. In (3_5), the first term on the right-hand side vanishes when {Z(e(k)] is noise-free; the second term lumps together all other errors present in {U(k)}. To simplify the presentation, let ?Zl-‘& = V be expressed as

for any appropriate from (35) that

matrices 1

g(z) - [Hodz)‘Q~(z) Q(z) =

3 and (e Then it foliows + %z)‘QWz)

+ ~(z)‘Q?~tz)l[&z) ii(z)'Qfi(z) R +

+ ~11 1

+ a(~)[ yf

+9

@Z)].

This in turn results in the following modification cIosed-loop stability criterion:

(38) of the

(39) where @(lr)llr_and n(z) and &(z)Z(z) (39) is computable. known.

(34)

When the data used are noise-free and no data in 0 are discarded, N = 0. Second, we assume that the data of (x((k)} are the corrupted data by n_oise-corrupted, and denote E(k) = x((k) + Z(k), where {5(k)} may be a transducer error or a disturbance signal: We account for corrupting noise by replacing E(‘(k) with x((k) in the computation of DII, Z and &$,from (15). Because of the presence of {Z(k)}, the data sequence computed from (15) will differ from the true fir(k). Let H(z)denote the z transform of this data sequence, and S(r) and T(z) the z transforms of (a(k)} and {T(k)] respectively. From e5), we infer that H(z) = H,,(z) + s(z), where z(z) = s(z)7(z))‘. In order to examine the error in {U(k)}, let the z transjorm of the comput_ed sequence of {U(k)} be denoted by U(z), and separate U(z) into U(z) = %(z) + G(z), where @U(z)is the true optimal control input and n(z) is the z transform of the error in {U(k)}. Then, if we denote the z transform of {N(k)} by N(z), the following z-domain equation of (34) is inferred: N(z) = M&)

k(z) = +(z)G(z))l

IIW(k)llF are the Frobenius norms of respectively, at z = e(21rrr)kl,Note that even though the plant model is not

4. lil~ir~ii~e example The example plant is defined by the following (A, B, C) triplet:

= HlfQf + (H’QW + iR)aj, where N(k) is m X m.

In general, the Des of IN(k)}, (F(k)) and {U(k)} serve as good estimates of N(z), H(z) and U(z) respectively for a11z on the_ unit circle. Furthermore, s(z) can be estimated from IT-‘, because T(z) is known. Hence a computation of n(z)/:=ae can be carried out, provided that an estimate of the noise spectrum, 3(z)/,=,.jo, is available. The presence of n(z) in D(z) also causes an erfor to appe.a_ in Y(z). From (19), we see that the error in Y(z) is D(z)Z’(z), where Z(z) is the z transform of {g(k)}. These errors induce additiona errors into (31), By using (29), we make the following mod~cation in (31):

t36J

0.7

r

0.5

0

0

rl

1

1

i

A=

c=[P

H.,

A,

:.,1.

In this exampie, these matrices are assumed not to be known explicitty. Design of the LQ regulator proceeds with the following steps. (i) Generate an open-loop system step response, sampling a total of 200 x(k). (ii) Compute {H(k)) from {Z’((k)}using (15), with T(k) = I. (iii) Compute {U(k)} from (14), using design parameters of p= IO-“, Q =I?= I and p= 100. It is found that I/fl(/c)ll

50, and therefore this is a steady-state solution. (iv) Compute {Y(k)} from {f./(k)}, using (19). Then the computation of 4 is carried out, using various values of Pand co~es~nding g(z) = z’in a trial process. For e = 1, the results are as follows: -0.71701 0.16342 0.44695 -7.22212 11.47523I ’ [ -12.01 p = 0.6071 0.19716 -0.4377 1 5.67504 1.6051 -3.56325 ’ [ I

F;, =

This is a noise-free soiution; hence E = 0. A test with noisey data is also conducted. In this test, unbiased and random white sensor noise is added to the data of {X(k)}. The root-mean-square noise-to-signal ratio of the data is 1%. It is found that the noise-corrupted {H(k)} decays down to the noise level for k > 70. Hence a value of k = 100

Brief Papers

Astrom, K. J. and B. Wittenmark (1989). Adaptive Control. Addison-Wesley, Reading, MA. Chan, J. T. (1995). Multivariable controller for redundantly actuated systems-A numerical design approach for open-loop data. Automatica, 35,781-786. _~_________________________________________ ‘$1(k) ,.--Chen, C. T. (1984). Linear System Theory and Design. Holt, Rinehart, and Winston, New York. Doyle, J. C. and G. Stein (1979). Robustness with observers. IEEE Trans. Autom. Control, AC-24,607-611. Franklin, G. F. and J. D. Powell (1980). Digital Control of Dynamic Systems. Addison-Wesley, Reading, MA. Golub, G. H. and C. F. Van Loan (1989). Matrix ‘w(k) Computations, 2nd ed. Johns Hopkins University Press, i..--.‘~ ________________ _--______________ ______________ b, ,(a) Baltimore. Ljung, L. and T. Siiderstriim (1983). Theory and Practice of ‘h(k) Recursive Identi$cation. MIT Press, Cambridge, MA. Ljung, L. (1987). System Identification-Theory for the User, Prentice-Hall, Englewood Cliffs, NJ. Patel, R. V. and N. Munro (1982). Multivariable System \.-~-~-.-._._._._ ‘b(k) Theory and Design. Pergamon Press, Oxford. Reid, J. G. (1983). Linear System Fundamentals, Continuous and Discrete, Classic and Modern. McGraw-Hill, New York. Sage, A. P. and C. C. White III (1977). Optimal System Control. Prentice-Hall, Englweood Cliffs, NJ. I I I I I 20.0 40.0 60.0 80.0 too.0 Appendix-Recursive computation of the optimal input Discrete time k N(k) =

\ \/

I.0

407

41,1(k) ‘h,,(k) h(k) b&k) h(k) b(k)

I

Fig. 1. Closed-loop system simulation. is selected. The design parameters of this test are the same as those of the noise-free test, but {U(k)} is computed by a recursive formulation in this case. Note that the initial value of 4 used is 9, = 0. A first-order controller results from this test, and shown as follows: -0.69936 F;, = [ -11.63871

EI =

sequence

If we let U(k) = 0 for all k > that

(mm where

+ W)cJ=

then it follows from (14)

-im@,

(A.1)

:I H(2) ‘+j= ; I1 [H(l)

1’ 1

0.17262 0.43114 -7.02872 11.14296

0.68902 0.19573 -0.43329 [ 5.50801 1.57878 -3.47598

k

H(l)

-”

A stability check of this solution using (39) results in the inequality

,

H(p’-k)

H(P + 1)

The closed-loop system step response of this controller is shown in Fig. 1 as six SISO data curves. The data matches the result produced by a model-based state feedback design. 5. Concluding remarks A data-based technique for linear-quadratic optimal regulator design has been presented. The method can be applied to MIMO optimal regulator design when the plant is open-loop stable, has more outputs than inputs and is left-invertible. The major advantages of the method are as follows.

Divide w into p block matrices, each block having dimensions of ti x pm. If we denote these block matrices by @, for i = 1, , p then (A.l) may be rewritten as f&=-Fp,$@~QH(i+l).

F,=(R+,z@:Q@,)m’J 64.2)

where U$ represents the estimate estimate of 111at time p + 1 will be

&+, = --&+I

(i) The desirable properties of the LQ design are retained.

of tiJ at time p. The

I$: (PTQH(i + 1).

(4.3)

(ii) The regulator is in an output feedback format, so that no state estimation is needed for implementation.

Equation (A.3) is a typical least-squares estimation formula. Its recursive form is given by

(iii) The computation of this optima1 output feedback controller does not require explicit identification of the plant model.

ojp+l = 4 - Fp+,@;+,g,,+i,

(iv) A reduced order controller, if one exists, may be synthesized during the design process without additional effort to reduce the plant model. References

Anderson, B. D. 0. and C. R. Johnson Jr (1982). Exponential convergence of adaptive identification and control problems. Automatica, 18,1-13.

(A.4)

where $ p+, = Q[H(p + 2) + @r+,~J~

*p = [H(P).

H(P -WI.

Finally, the following recursion formula for Fp is obtained from the matrix inversion lemma: F ,,+, = Fp - F,@;+I(I

+

Q~~+,F,~~+,)-‘Q~~+,F,, (A.3

with the initial condition Fo = R-‘.