Quasi-direct pole placement for time delay systems applied to a heat transfer set-up Tom´ aˇ s Vyhl´ıdal ∗ Wim Michiels ∗∗ Pavel Z´ıtek ∗ ∗ Center for Applied Cybernetics, Dept. of Instrumentation and Control Eng., Faculty of Mechanical Engineering, Czech Technical University in Prague, Technick´ a 4, 166 07 Praha 6, Czech Republic, (e-mail:
[email protected],
[email protected]) ∗∗ Department of Computer Science, Katholieke Universiteit Leuven, Celestijnenlaan 200A, B-3001 Heverlee, Belgium, (email:
[email protected])
Abstract: A method for quasi-direct pole placement is presented for a class of retarded time delay systems with a single control input. The pole placement is done in two steps. First, a number of dominant poles, smaller than the number of controller parameters, are directly assigned. This leads to constraints on the controller parameters. By using the singular value decomposition a simple parameterization of all controllers satisfying these constraints is obtained. In the second step, the remaining degrees of freedom in the parameter space are used to shift the remaining part of the system spectrum as far to the left as possible. It is done using an optimization procedure for nonsmooth, nonconvex functions. An extensive real-plant example is presented to demonstrate the application of the novel pole placement procedure.Copyright c IFAC 2009 ° Keywords: Retarded system, pole placement, singular value decomposition, nonsmooth optimization 1. INTRODUCTION The conventional state feedback scheme offers a unique opportunity to design fully the dynamics of a linear delay free system. The well known method for deriving the feedback gain is the pole placement method. For n − th order SISO systems, all n poles are assigned into desired positions and consequently n feedback gain coefficients are computed as a single solution in the parameter space, see e.g. Kautsky and Nichols, (1985). As it has been shown in Z´ıtek, (1997), Z´ıtek and Vyhl´ıdal, (2000), Michiels, et al., (2002), Vanbiervliet, et al., (2008), the classical state feedback can also be applied to adjust the dynamics of time delay systems. However, it has considerable limitations induced mainly by the infinite system spectrum and the limited degree of freedom in the parameter space. We consider a retarded system Z T Z T x(t) ˙ = dA(τ )x(t − τ ) + dB(τ )u(t − τ ) (1) 0
0
where x ∈ Rn is the vector of state variables, u ∈ Rp is the vector of system inputs, τ is the delay variable, which is bounded by 0 ≤ τ ≤ T . The functional matrices A(τ ), B(τ ) of compatible dimension describe the distribution of the delay; their discontinuities correspond to the system lumped delays. Let us mention that the Stieltjes ? This research has been supported by the Ministry of Education of the Czech Republic under the Project 1M0567, and by the Programme of Interuniversity Attraction Poles of the Belgian Federal Science Policy Office, under the project IAP P6-DYSCO and by OPTEC, the Optimization in Engineering Center of the K.U.Leuven.
integral form in (1) covers both the lumped and distributed delays, Gorecki, et al. (1989), Z´ıtek, (1997). The state feedback controller is considered in the classical form uc (t) = −K T x(t) (2) where uc is the single control input and K = [k1 , k2 , ...., kn ]T is the feedback gain. The remaining p − 1 inputs of the system are considered as the system disturbances. The stability of feedback system (1)-(2) is determined by the solutions of the characteristic equation M (λ, K) = det(sI − A(λ) + B(λ)K T ) = 0 (3) where Z T
A(λ) =
exp(−λτ )dA(τ ) Z 0T
B(λ) =
exp(−λτ )dB(τ ). 0
System (1)-(2) is exponentially stable iff all the infinitely many roots of equation (3) are located in the open left half plane, Hale and Verduyn Lunel, (1993), Michiels and Niculescu, (2007). As the time delay system is retarded, the following spectrum distribution features hold. Property 1. For any β ∈ R, the number of roots of (3) with <(λi ) > β is finite. Property 2. For |λi | > c >> 0, c ∈ R+ the roots of a retarded system tend to match the finite number of exponential curves with a positive exponents, Bellman and Cooke, (1963), Vyhl´ıdal and Z´ıtek, (2009) . A consequence of Property 1 is that a retarded system may have finitely many unstable roots only. A consequence of
the Property 2 is that infinite chains of roots converge to −∞ in real and ±∞ in imaginary part. As a consequence of both the properties, the dynamics of a retarded system are determined by a small number of rightmost roots, Z´ıtek, (1997). As it has been shown in Z´ıtek, (1997), Z´ıtek and Vyhl´ıdal, (2000), n system poles can be placed to desired positions by n feedback parameters in the same way as in the case of a delay free system. However, such a placement is effective if and only if the rest of the system spectrum is located to the left of these assigned dominant poles, which is not guaranteed. In fact, the direct pole placement relies on rather heuristic trial and error placement of the dominant poles. If the assigned poles are not isolated from the rest of the spectrum, the pole placement procedure is repeated with a different selection of the poles and so on. A fully automated pole placement based control design - continuous pole placement has been introduced in Michiels, et al., (2002) for retarded and in Michiels and Vyhl´ıdal, (2005) for neutral systems. The aim of this method is to shift the rightmost roots of the feedback system as far to the left as possible via small-step-wise numerically based shifting of the rightmost roots. In this way, the spectral abscissa α = sup(<(λi )), i = 1, 2, .., ∞, is minimised. More recently an optimization approach has directly been used to minimize the spectral abscissa α in Vanbiervliet, et al., (2008). As the spectral abscissa is a nonsmooth function, the gradient sampling algorithm was used, see Burke, et al., (2005), see also Burke, et al., (2006). Even though the minimization of the spectral abscissa α is very useful from the point of system stabilization, it is not a true pole placement method in fact. The resulting system dynamics is determined by a number of more than n poles located close to min α. Although the dynamics determined by these poles is the fastest possible, there is no guarantee of proper damping of the system modes as the imaginary parts of the poles are not assigned. In this paper, we provide a pole placement method which combines the positive features of direct and optimization based pole placement. The aim is to place less than n system poles and to push the rest of the spectrum as far to the left as possible via utilizing the remaining degree of freedom in the parameter space. In this way, we want to guarantee that the assigned poles will be the dominant poles of the feedback system determining the system dynamics. 2. STABILIZATION WITH POLE LOCATION CONSTRAINTS Due to the linearity of the characteristic function of system (1) with respect to the parameters of K, it can be transformed into the following form n X M (λ, K) = Q(λ) + ki Pi (λ) = 0 (4) k=1
where Q(λ) = det(λI − A(λ)) and Pi (λ) = 1..n are quasi-polynomials.
∂M (λ,K) ,i ∂ki
=
Assigning a real pole to c yields the following constraint on the gain values:
n X
Pi (c)ki = −Q(c).
k=1
Similarly, assigning a complex conjugate pair of poles, c ± di, results in n X <(Pi (c + di))ki = −<(Q(c + di)), k=1
n X =(Pi (c + di))ki = −=(Q(c + di)). k=1
In this way, assigning 0 ≤ m ≤ n poles to λ1 , . . . , λm eventually results in set of m constraints which can be written in the form SK = R, (5) where S ∈ Rm×n and R ∈ Rm×1 . If the pole assignment problem is solvable, that is, equation (5) has solutions, then the latter describes a manifold in the parameter space corresponding to the presence of λ1 , . . . , λm in the spectrum of feedback system (1)-(2). Let us simplify this manifold using another basis of the parameter space. Using the singular value decomposition S = U ΣV ∗ , where (·)∗ denotes the complex conjugate transpose, condition (5) becomes ΣL = U ∗ R (6) where L = V ∗ K. (7) Letting L = [l1 , . . . , ln ]T and U ∗ R = [¯ r1 , . . . , r¯m ]T , Σ = T [σ1 , ...σm ] , and assuming that S is of full (row) rank, we finally get the following expression for computing the parameters of li , i = 1..m corresponding to the assigned poles λi , i = 1..m: l1 = r¯1 /σ1 l2 = r¯2 /σ2 . (8) .. . lm = r¯m /σm The following can be concluded: Proposition 3. Consider the feedback system (1)-(2). By assigning m system poles λ1 , . . . , λm to its characteristic equation (4), the system of m equations (5) is obtained, constraining the feedback gain K. If S has full row rank, then the constraint (5) is equivalent to (8), where the relation between K and L = [l1 , . . . , ln ]T is described by L = V ∗ K, K = V L. (9) Proposition 4. Controller parameterization. Assume that S has full row rank. Considering parameters l1 , . . . , lm as fixed according to (8), and lm+1 , . . . , ln being kept as free parameters, one obtains a parametrization of all controllers that assign m poles to λ1 , . . . , λm . The next step of the pole placement procedure is given as follows. Consider m poles are assigned as described in Propositions 3-4. Thus, the parameters l1 , . . . , lm are fixed and the parameters lm+1 , . . . , ln are available for further adjustment of the system spectrum. Applying the algorithm described in Vanbiervliet, et al., (2008), parameters lm+1 , . . . , ln are to be used to push the other rightmost poles as far as possible to the left. However, unlike in Vanbiervliet, et al., (2008), where the spectral
abscissa α = sup(<(λi )), i = 1, 2, .., ∞, is minimised, we consider to minimize the function α ¯ (lm+1 P ¾ ½ , . . . , ln ) = Q(λ) + ki (lm+1 , . . . , ln )Pi (λ) = 0 . sup <(λ) : Πm i=1 (λ − λi ) (10) Technically, the evaluation of α ¯ is quite straightforward. First, rightmost zeros of M (λ, K) given in (4) are computed either directly by the rootfinder recently presented in Vyhl´ıdal and Z´ıtek, (2009) or as the rightmost eigenvalues of system (1)-(2) by e.g. DDE-BIFTOOL of Engelborghs, et al. (2002). In the next step, the (invariant) eigenvalues λ1 , . . . , λm are removed from the spectrum and the objective function is obtained as α ¯ = supi>m (<(λi )). Analogous as in Vanbiervliet, et al., (2008), it can be shown that that the function α ¯ is smooth almost everywhere (i.e., except for a set in the parameter space with measure zero). More precisely, the function is smooth whenever there is only one non-assigned eigenvalue whose real part is equal to α ¯ . In such a case one can compute the gradient of α ¯ from the sensitivity of this eigenvalues w.r.t. the parameters. The gradient sampling algorithm, presented in Burke, et al., (2005) and used in Vanbiervliet, et al., (2008) can now be applied to minimize α ¯ (lm+1 , . . . , ln ). The algorithm, developed for the minimization of non-smooth, non-convex functions, which are differentiable almost everywhere, uses a technique referred to as gradient sampling to construct a gradient bundle. The algorithm is essentially the same as the Steepest Descent method, apart from the fact that a non-smooth steepest descent direction is computed. This can be defined as − arg min kzk (11) ¯c) z∈∂c φ(L ¡ ¢ ¯ c denotes the generalised gradient (referred where ∂c φ L ¯ c , where L ¯c = to as the Clarke Subdifferential ) at L [lm+1 , . . . , ln ]. It is given by ½ ¾ ¡ ¢ ¡ ¢ ¯ ¯ ∂c φ Lc = conv lim ∇¯ α L , (12) ¯ L ¯c L→
where conv(.) denotes the convex hull function. Note that the generalized gradient reduces to the classical notion of gradient whenever the function is differentiable. An application of the gradient sampling algorithm leads to a monotonic decrease of the objective function until a local minimum is¡ reached. Such a minimum is characterized ¢ ¯ c . In order to speed up the optimization by 0 ∈ ∂c φ L process one can start by applying the BFGS algorithm, which performs surprisingly well on nonsmooth problems (see Lewis and Overton, (2009) for a deep analysis), and switch to the gradient sampling algorithm when the BFGS algorithm stagnates. The pole placement procedure can now be summarised as follows Algorithm 5. Quasi-direct pole placement. Consider feedback system (1)-(2), the pole placement is done in the following steps (1) Select poles λ1 , . . . , λm , m < n to be assigned. (2) Compute the parameters l1 , . . . , lm as described in Propositions 3-4. (3) Minimize the spectral abscissa α ¯ over the parameters lm+1 , . . . , ln .
(4) If min α ¯ < min <(λi ), i = 1..m, accept the result and trasform L to K by (9), otherwise select different λ1 , . . . , λm , m < n and go to step 1. Remark 6. According to Step 4 of Algorithm 5, after assigning the poles and running the optimization procedure, one has to check whether the assigned poles are located to the right of the current α ¯ . If not, the dominance of the assigned poles is not guaranteed and the pole placement procedure must be run again with a different setting. Remark 7. In order to estimate a reasonable distance from the imaginary axis for assigning the poles, the spectral abscissa α(K) can be minimised over all parameters available in K in the first place, as it is described in Vanbiervliet, et al., (2008). It is then advised to assign poles such that min α(K) < <(λi ) < 0, i = 1, .., m. (13) Indeed, note that if the rightmost assigned pole has real part smaller than min α(K), then the assigned poles cannot be the rightmost ones. Remark 8. A good starting value for the gradient sampling algorithm can be obtained as arg min {kK − K0 k : SK = R} , (14) where K0 = arg min α(K). In words, the expression (14) corresponds to the gain value in the manifold (5) closest to K0 . In this way, it is expected that the local minimizer of α ¯ , found by Algorithm 5, also corresponds to a gain value close to K0 , which is important to prevent a large qualitative change in the uncontrolled part of the spectrum. An explicit computation of (14) leads to lm+1 l1 .. .. † . = [vm+1 · · · vn ] K0 − [v1 · · · vm ] . , ln lm (15) where V = [v1 · · · vm ], l1 , . . . , lm are given by (8), and (·)† denotes the Moore-Penrose inverse. 3. POLE PLACEMENT CONTROL DESIGN OF EXPERIMENTAL HEAT TRANSFER SET-UP The above control design has been tested on a model of a real plant - experimental heat transfer set-up. The set-up, which is shown in Fig. 1 and which has been described comprehensively in Vyhl´ıdal and Z´ıtek, (2009) consists of two closed heating circuits with forced circulation of water as the heat transfer medium. Both circuits are equipped with heathers and coolers. The heat exchange between the circuits takes place in a heat exchanger. These main subsystems of the set-up are connected by pipe-lines. Obviously, the rations between the lengths of pipe-lines and flow velocities determine the main transport delays in the heat transfer throughout the system. Besides, the set-up is equipped with three mixing valves for controlling the flows in some parts of the set-up. The control inputs of the system are the performances of heaters and coolers and the system outputs are temperatures measured at inputs and outputs of the main subsystems. In Vyhl´ıdal and Z´ıtek, (2009), the linear model of the setup was proposed, which consists of ten linear differential equation in the form of (1), with the following L-transform of the system matrices
Table 1. Results of the quasi-direct pole placement. SN - Setting Number, λi - assigned poles, min α ¯ - minimum of the spectral abscissa of the rest of the spectrum, ki , i = 1..11 coefficients of the feedback gain in (17). SN λi
1 -
2
3
4
-0.01
−0.03 ± 0.03i
-0.025, -0.035 −0.03 ± 0.03i
min α ¯
-0.0565
-0.0632
-0.0576
-0.0724
k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11
-5.4349 3.5879 -1.4411 -3.7043 24.6167 -2.1778 9.6924 -4.5121 -14.6312 11.3514 -0.7562
-0.0579 8.5947 -1.6531 -3.8707 32.9189 3.5683 10.5949 -6.1505 -21.5687 4.0466 -0.2690
-3.5278 5.8479 -2.0390 -3.1424 26.2714 -1.3073 9.4212 -4.7706 -18.1263 7.1200 -0.4747
-0.8946 8.6153 -6.0585 -4.5942 33.7005 -0.8714 6.9169 -6.0465 -23.5958 2.9531 -0.1812
form
uc (t) = −K T x ¯(t) (17) where x ¯(t)) = [x(t) I(t)]T is the extended state variable vector where Zt I(t) = (ySET (η) − y(η))dη (18) T
Fig. 1. Experimental heat transfer set-up with heat exchangers. a1,1 (λ) ... a1,11 (λ) .. .. .. , A(λ) = . . . a11,1 (λ) ... a11,11 (λ) (16) b1,1 (λ) b1,2 (λ) .. .. B(λ) = . . b11,1 (λ) b11,2 (λ) , a2,1 (λ) = where a1,1 (λ) = − 51 , a1,6 (λ) = 0.973exp(−λ15) 5 0.983exp(−λ5) 0.96exp(−λ23) 1 , a2,2 (λ) = − 25 , a3,2 (λ) = , 25 5 1.425 , a (λ) , a (λ) a3,3 (λ) = − 51 , a4,3 (λ) = 0.575 = − = 4,4 4,7 3 3 0.975exp(−λ3) 0.85 0.85 , a5,5 (λ) = 23 , a4,8 (λ) = 23 , a5,4 (λ) = 5 0.9exp(−λ5) 1 1 , a6,6 (λ) = − 17 , a7,7 (λ) = − 5 , a6,5 (λ) = 17 exp(−λ29) 0.85 a8,3 (λ) = 23 , a8,4 (λ) = 0.85 − 63 23 , a8,7 (λ) = 0.97exp(−λ5) 1.425 0.575 , a (λ) , a (λ) , a9,9 (λ) = = − = 8,8 9,8 3 3 5 0.92exp(−λ5) 1 1 , a10,10 (λ) = − 15 , b7,1 (λ) = − 5 , a10,9 (λ) = 15 0.042exp(−λ6) exp(−λ7) , b10,2 (λ) = , while the other terms in 63 15 the matrices A(λ) and B(λ) are equal to zero. The system states are the temperatures measured at different places of the set-up. We consider here only two out of four system inputs from Vyhl´ıdal and Z´ıtek, (2009) u = [uc , ud ]T . The first - control signal uc is the setpoint value of the slave PI control loop that controls the performance of the left heater (the slave control loop is included in the model of the set-up). The second signal ud adjusts the performance of the left cooler and it is considered here as the disturbance variable. The controlled variable y is the temperature measured at the output of the left cooler, which is the last state variable in the vector x. The feedback controller is considered in the
0
is an additional state variable introduced for the control purposes, ySET being the set-point of y. In Fig. 3, the rightmost poles of the open loop system are shown. As can be seen, due to the astatism introduced by (18) one pole is placed at the complex plane origin. The remaining poles determine the modes of the set-up model. The system step responses to both system inputs can be seen in Fig. 2. The parameters of the feedback gain has been designed by the quasi-direct pole placement described in Algorithm 5. The results of the synthesis are given in Table 3 where the assigned poles are displayed as λi in the first row. In the second row, the spectral abscissa of the remaining part of the spectrum is given and the resulting coefficients are given below. As can be seen, no poles are assigned for the setting SN 1, which presents the results of minimizing the spectral abscissa α over all 11 parameters of the feedback gain K using the algorithm proposed in Vanbiervliet, et al., (2008). The resulting spectrum for this setting is shown in Fig. 4 and the corresponding set-point and disturbance rejection responses are shown in Fig. 2. In fact, this setting corresponds to the fastest possible response by (17) as the poles are shifted as far to the left as possible. As it has been mentioned in Remark 7, it is suitable to place the poles to the right of min α = −0.0565 of SN 1 in Table 3. First, only a single pole has been assigned λ1 = −0.01 resulting in min α ¯ = −0.0632 (SN 2 in Table 3, see the spectrum in Fig. 5. As can be seen in Fig. 2, the corresponding responses are fairly slow. Secondly a pair of poles λ1,2 = −0.03 ± 0.03i has been assigned resulting in min α ¯ = −0.0576 (SN 3 in Table 3, see the spectrum in Fig. 6 ) and fast and well dumped responses,
0.25
4.5
3
1 3
4
2
3
4
3.5
0.2 0.15
1
0.1
4
0.05
1
0
3 2.5
ℑ(λ)
−0.06 −0.03
y
2
2
2
Step
1.5
−1
ySET
1
−2
0.5 0
0
0
−3 −0.2
Step
−0.15
−0.1
−0.05
0
ℜ(λ) 0
200
400
600
800
1000
1200
t
Fig. 2. Step responses of the set-up model (Step), uc = 4 at t = 100 and ud = 0.5 at t = 700 and the setpoint and disturbance rejection responses of the feedback system with settings (SN) 1..4 from Table 3, ySET = 3 at t = 100 and ud = 0.5 at t = 700.
Fig. 5. Spectrum of the feedback system for SN 2 in Table 3. 0.25
3
0.2 0.15
2
0.1 0.05
1
0.25
3
ℑ(λ)
0.15
2
0
0
0.1
−1
0.05
1
0 −0.06 −0.03
ℑ(λ)
0 −0.06 −0.03
0.2
0
−2
0 −3 −0.2
−1
−0.1
−0.05
0
ℜ(λ)
−2 −3 −0.2
−0.15
Fig. 6. Spectrum of the feedback system for SN 3 in Table 3. −0.15
−0.1
−0.05
0
0.25
ℜ(λ)
3
Fig. 3. Spectrum of the open loop system.
2 1
0.2
0 −0.06 −0.03
ℑ(λ)
0.15
2
0.1 0.05
0.25
3
0.2 0.15
0.1
0
0
0.05
1
−1
0
ℑ(λ)
−0.06 −0.03
0
−2
0
−3 −0.2
−1
−0.15
−0.1
−0.05
0
ℜ(λ) −2 −3 −0.2
−0.15
−0.1
−0.05
0
ℜ(λ)
Fig. 4. Spectrum of the feedback system for SN 1 from Table 3. however with a considerable overshoot. Finally, four poles have been assigned given in SN 4 of Table 3 resulting in min α ¯ = −0.0724, see see the spectrum in Fig. 7. As can be seen in Fig. 2, the responses of SN 4 are the best achieved - fast, well dumped and with small overshoot.
Fig. 7. Spectrum of the feedback system for SN 4 in Table 3. Let us describe control synthesis e.g. for SN 4 in more detail where m = 4 have been assigned. First, the singular value decomposition is applied over the system of equation (5). Then, the parameters li , i = 1..4 are obtained directly from (8) when assigning the four selected poles. The remaining parameters are obtained as the results of optimization approach of minimizing α ¯ given in (10). The results of the optimization procedure over the parameters are shown in Fig. 8 starting from li = 0, i = 5..11. The
REFERENCES
α (l5,..,l11)
−0.02 −0.04 −0.06 −0.08
0
50
100
150 200 250 iteration number
300
350
300
350
30 l5
20
l7
l10
5
l ,..,l
11
10
l6
0 −10 l9
−20 −30
0
50
100
l11
150 200 250 iteration number
l8
Fig. 8. Evolution of α ¯ and the parameters li , i = 5..11 as the result of the optimization process for locating min α ¯ in the parameter space. The initial value of the parameters is obtained from (15). The setting is SN 4. optimization procedure is stopped at the iteration 465 when the min α ¯ is reached. Finally, the parameters of the gain K are obtained by (9).
4. CONCLUSION A novel pole placement procedure for retarded time delay systems with a single control input has been proposed. It combines direct and optimization based pole placement. The aim is to place a small number of distinct dominant poles determining the dynamics of the system. These poles should be placed sufficiently far from the imaginary axis to guarantee both the robust stability of the system and sufficiently fast dynamics. Besides, the ratio between real and imaginary parts of the poles should be large enough in order to guarantee a proper damping in the system responses. Further on, an optimization based pole shifting algorithm is used to push the remaining rightmost poles as far to the left as possible. Only such pole distribution is acceptable in which the assigned poles are sufficiently isolated from the rest of the spectrum. In order to demonstrate the capability of the approach, the pole placement has been applied to a control synthesis of a real plant model - experimental heat transfer set-up. Besides the demonstration of the novel pole placement procedure, the contribution of the example section lies in a new application of optimization based fixed structure control design techniques. Finally, although the control design technique was illustrated by means of static instantaneous state feedback the results can be trivially extended to the case where the control law includes multiple delayed measurements or even distributed delays.
Engelborghs K., T. Luzyanina, and D. Roose, (2002) Numerical bifurcation analysis of delay differential equations using DDEBIFTOOL, ACM Trans. Math. Softw., 28, pp. 121. Bellman, R., and Cooke, K.L. (1963), Differentialdifference equation, Academic Press, New York. Burke, J., A. Lewis, and M. Overton, (2005), A robust gradient sampling algorithm for nonsmooth, nonconvex optimization, SIAM J. Opt., 24, pp. 567584. Burke, J., D. Henrion, A. Lewis, and M. Overton, (2006) HIFOO - A matlab Package for Fixed-Order Controller Design and H-infinity optimization, in Proceedings of ROCOND 2006, Toulouse, France. Gorecki, H. Fuksa, S., Grabowski, P. and Korytowski, A., (1989),Analysis and Synthesis of Time-Delay Systems. Polish Scientific Publications, Warsaw Hale, J. K., and Verduyn Lunel, S. M., (1993), Introduction to functional differential equations, vol. 99 of Applied Mathematical Sciences, Springer Verlag New York Inc. Kautsky, J. and N.K. Nichols, (1985), Robust Pole Assignment in Linear State Feedback, Int. J. Control, 41, pp. 1129-1155. Michiels, W. and Niculescu, S., (2007), Stability and Stabilization of Time Delay Systems: An Eigenvalue Based Approach. Advances in Design and Control, SIAM. Michiels, W., Engelborghs, K., Vansevenant, P. and Roose, D. (2002) Continuous pole placement method for delay equations, Automatica 38(5), pp. 747-761. Michiels, W. and Vyhl´ıdal, T, (2005), An Eigenvalue Based Approach for the Stabilization of Linear Time-Delay Systems of Neutral Type, Automatica, Vol. 41, pp. 991998 Lewis, A.S. and Overton, M.L., (2009), Nonsmooth Optimization via BFGS, SIAM Journal on Optimization. Submitted. Vanbiervliet J., K.Verheyden, W. Michiels, S.Vandewalle (2008), A nonsmooth optimization approach for the stabilization of time-delay systems ESAIM Control, Optimisation and Calcalus of Variations, Vol.14, No.3, pp.478-493. Z´ıtek, P., (1997) Frequency Domain Synthesis of Hereditary Control Systems via Anisochronic State Space. Int. Journal of Control, Vol. 66, No. 4, pp. 539-556. Z´ıtek, P. and Vyhl´ıdal, T., (2000) State feedback control of time delay system: conformal mapping aided design. In: 2nd IFAC Workshop on Linear Time Delay Systems. Ancona, Italy, pp. 146-151. Vyhl´ıdal, T. and Z´ıtek, P. (2009), Mapping based algorithm for large-scale computation of quasi-polynomial zeros, IEEE Transactions on Automatic Control, Vol. 54, No. 1, pp. 171-177. Vyhl´ıdal, T., Z´ıtek, P., Paul˚ u, K. (2009) Design, modelling and control of the experimental heat transfer set-up, In Topics in Time Delay Systems Analysis, Algorithms, and Control, edited by: R. Sipahi, W. Michiels. J. J. Loiseau, S.-I. Niculescu, Lecture Notes in Control and Information Sciences, Springer. To appear.