Dan Ma
email: madan@mail.neu.cn.edu

Islam Boussaada
email: islam.boussaada@l2s.centralesupelec.fr

Jianqi Chen
email: jianqchen2-c@my.cityu.edu.hk

Catherine Bonnet
email: catherine.bonnet@centralesupelec.fr

Silviu-Iulian Niculescu
email: silviu.niculescu@l2s.centralesupelec.fr

Jie Chen
email: jichen@cityu.edu.hk

PID Control Design for First-Order Delay Systems via MID Pole Placement: Performance vs. Robustness

Keywords: Time-delay systems, PID control, delay robustness, pole placement, steady-state tracking

The PID control is favored in controlling industrial processes for its ease of implementation. In this paper we present an analytical study, with an aim at stability robustness, tracking performance, and pole placement of first-order unstable plants under PID control. We employ a multiplicity-induced dominancy (MID) strategy to design and tune PID controller gains, which is shown to possess desirable properties that ensure specified decay rate of the closed-loop response and the stability robustness of the closed-loop system in spite of variations in the delay parameter. We also study the tradeoff between delay robustness and tracking performance of PID controllers. We show that under the constraint of steady-state tracking, the maximal delay robustness range can be computed by solving a unimodal pseudo-concave optimization problem.

INTRODUCTION

A canonical case of study for PID (proportional-integralderivative) control is the first-order delay system

P τ (s) = K 1 + T s e -τ s ,
where K > 0 is the plant gain, τ ≥ 0 represents the time delay in the controlled process, and T ≥ 0 is referred to as the time constant of the system. Traditionally, the implementation of a PID controller is done via parameter tunning, based on such empirical rules as Ziegler-Nichols method and its offsprings [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF]. The PID control has been long

well-known for its simplicity, ease of implementation and cost-effectiveness. Even in a post-modern and information-centric era of present time, it remains to stand out as a most favored control technology and dominates industrial control systems; by documented statistics, more than 95% of the physical-layer control loops are still enabled by PID controllers, with proven, widespread acceptance by the industrial control community [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF][START_REF] Samad | A survey on industry impact and challenges thereof[END_REF].

For its recognized advantages, there has been sustained interest in PID control. Some of the recent work on PID control has been focused on analytical studies, which have led to new findings, improved design and tuning rules on e.g., time delay systems, nonlinear systems, and extensions to multi-agent systems [START_REF] Xu | Input delay margin for consensusability of multi-agent systems[END_REF][START_REF] Ma | Delay consensus margin of first-order multiagent systems with undirected graphs and PD protocols[END_REF]. Specifically, steady-state tracking performance was examined in [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF]; [START_REF] Silva | New results on the synthesis of PID controllers[END_REF], the delay robustness problems were studied in [START_REF] Silva | New results on the synthesis of PID controllers[END_REF]; [START_REF] Ma | Delay margin of low-order systems achievable by PID controllers[END_REF]; [START_REF] Chen | Exact computation of delay margin by PID control: It suffices to solve a unimodal problem![END_REF], and nonlinear regulation problems were considered in [START_REF] Zhao | PID controller design for second order nonlinear uncertain systems[END_REF]. Of particular relevance to this paper, there is a long stream of works concerning design problems such as pole placement of timedelay systems by PID controllers, based on tools ranging from algebraic techniques [START_REF] Wang | New results on PID controller design of discrete-time systems via pole placement[END_REF][START_REF] Emre | Digital PI-PD controller design for arbitrary order systems: Dominant pole placement approach[END_REF][START_REF] Das | Delay handling method in dominant pole placement based pid controller design[END_REF], numerical computation [START_REF] Michiels | Continuous pole placement for delay equations[END_REF][START_REF] Pavel | Dimensional analysis approach to dominant three-pole placement in delayed PID control loops[END_REF][START_REF] Firouzbahrami | Reliable computation of PID gain space for general secondorder time-delay systems[END_REF][START_REF] Jaromír | PID controller tuning via dominant pole placement in comparison with Ziegler-Nichols tuning[END_REF][START_REF] Zhang | An improved analytical tuning rule of a robust PID controller for integrating systems with time delay based on the multiple dominant pole-placement method[END_REF], to performance and robustness considerations [START_REF] Srivastava | A PI/PID controller for time delay systems with desired closed loop time response and guaranteed gain and phase margins[END_REF][START_REF] Tavakoli | Analytical PID control design in time domain with performancerobustness trade-off[END_REF][START_REF] Rao | Tuning optimal PID controllers for open loop unstable first order plus time delay systems by minimizing ITAE cri-terion[END_REF].

Also of close pertinence are studies on open-loop unstable systems with time delays. Physical processes that can be suitably modeled as an unstable delay plant are widely found in the process industry, including, e.g., bioreactors and tank fermenters (see, e.g, O'Dwyer ( 2013) and the references therein). For such systems, the delay margin [START_REF] Middleton | On the achievable delay margin using LTI control for unstable plants[END_REF] furnishes an intrinsic measure of robustness with respect to variations in time delay, defined as the maximal range of time delay within which the system can be stabilized robustly by a fixed controller. This notion has been extensively studied recently, based on different methods and techniques (see, e.g., [START_REF] Bresch-Pietri | Adaptive control scheme for uncertain time-delay systems[END_REF]; [START_REF] Gaudette | Adaptive stabilization of a class of time-varying systems with an uncertain delay[END_REF]; [START_REF] Foias | Robust Control of Infinite Dimensional Systems: Frequency Domain Methods[END_REF]; [START_REF] Ju | Further results on the achievable delay margin using LTI control[END_REF]; [START_REF] Miller | Stabilization in the presence of an uncertain arbitrarily large delay[END_REF]; [START_REF] Zhong | Duality System in Applied Mechanics and Optimal Control[END_REF]; [START_REF] Zhou | Properties of the parametric lyapunov equation-based low-gain design with applications in stabilization of time-delay systems[END_REF][START_REF] Zhou | Truncated predictor feedback for linear systems with long timevarying input delays[END_REF]). Analytical and computational results have been previously obtained for determining the delay margin achievable by PID controllers [START_REF] Silva | New results on the synthesis of PID controllers[END_REF][START_REF] Ma | Delay margin of low-order systems achievable by PID controllers[END_REF][START_REF] Chen | Exact computation of delay margin by PID control: It suffices to solve a unimodal problem![END_REF]. This fundamental limit, however, is attained with no regard to other control objectives, and hence is in inherent conflict with a system's transient and steady-state performance. In a more realistic design, one must then incorporate performance objectives and address the tradeoff between performance and robustness. This motivates our work in this paper.

We study the analytical design of PID controllers for delay robustness and performance. We consider first-order unstable plants. This consideration partly stems from the fact that industrial processes are often modeled by first-and generally low-order systems, partly due to the limitation of PID controllers in controlling high-order dynamics. While limited in scope, the study of first-order plants does allow an in-depth analysis, availing analytical results that render useful insights. We develop a poleplacement design via the so-called multiplicity-induced dominancy (MID), which seeks to place the dominant poles of the closed-loop system at a same location, resulting in a multiple dominant pole. The design consequently ensures such transient property as the degree of stability, or alternatively the decay rate of the closed-loop response. Of equal importance, it also renders a certain amount of delay robustness, with a guaranteed range of delay within which the plant can be robustly stabilized. Interestingly, the result also demonstrates a clear tradeoff between the delay robustness range and the decay rate, and in the limit the range approaches to the delay margin, i.e., the largest range possible to maintain robust stability, whilst the closed-loop system becomes marginally stable.

The analytical MID pole placement design for time-delay systems by itself poses a nontrivial problem, even for a first-order plant. The technical intricacy is further compounded by additional requirement to ensure delay robustness. Previously an analytical proof of the dominancy of a spectral value for a scalar state-space equation with a single delay was presented in [START_REF] Hayes | Roots of the transcendental equation associated with a certain difference-differential equation[END_REF].

The dominancy property was also established for scalar state-space delay equations in Boussaada et al. (2016), and further exploited in [START_REF] Boussaada | Further remarks on the effect of multiple spectral values on the dynamics of time-delay systems. application to the control of a mechanical system[END_REF][START_REF] Boussaada | Multiplicity-induceddominancy in parametric second-order delay differential equations: Analysis and application in control design[END_REF] for second-order retarded differential equations. For more general pole placement problems via PID control, numerical optimization-based methods developed in ( [START_REF] Vanbiervliet | A nonsmooth optimisation approach for the stabilisation of time-delay systems[END_REF][START_REF] Ramirez | Design of proportional-integral-retarded (PIR) controllers for second-order LTI systems[END_REF] incorporated as well the multiplicity-induced dominancy. It is useful to note that in all the aforementioned results, the assumption of a retarded system is critical. In an important distinction, however, with the use of PID control, the system under consideration in the present paper constitutes invariably a neutral delay system. Unlike retarded systems which admit a finite number of roots in any right half plane, neutral delay systems may have infinitely many roots located in the open right half plane. This makes our design considerably more challenging and it requires a nontrivial extension of the multiplicityinduced dominancy analysis to neutral delay systems.

We also investigate the tradeoff between the delay margin and the steady-state tracking performance. While in the limit the delay margin can be analytically determined for first-order systems, under performance constraint it generally poses a constrained parametric optimization problem which in turn requires a brute-force search in the three design parameters of the PID controller. Nevertheless, we show that under a prescribed specification of tracking error, the problem can be solved computationally as a unimodal pseudo-concave program, i.e., a univariate smooth convex optimization problem. As such, under the steady-state tracking objective, the maximal delay robustness range can be computed using convex programming methods, or gradient-based bisection methods. The attainable range is found to be inversely proportional to the tracking accuracy. The solution is made possible by recasting the problem as one of nonlinear programming, and tackled by employing the well-known Fritz John condition.

The remaining paper is organized as follows. In Section 2 we introduce the mathematical background required in the sequel, together with the notion and facts of delay margin. We then provide the multiplicity-induced dominancy design in Section 3, where we analyze the rightmost root of the closed-loop quasipolynomial in details.

The multiplicity and dominance of the rightmost root is established. In Section 4, we further establish the delay robustness of this multiplicity-induced dominancy design, showing that a range of delay exists that guarantees the closed-loop stability. Section 5 studies the tradeoff between tracking performance and delay robustness. We show that under the constraint of steady-state tracking, the maximal delay robustness range can be computed by solving a unimodal pseudo-concave optimization problem. An illustrative example is then provided in Section 6, and the paper concludes in Section 7.

Partial results of this paper were previously presented in Chen and Chen (2020); Ma et al. (2020a); the current paper includes updated results with full proofs.

2 Preliminary background

Mathematical preliminaries

The primary technical machineries to be used in the sequel dwell on complex variables theory and mathematical programming problems. Consider the complex function

G(s) = r(s) Q 0 (s) + Q τ (s)e -sτ , (1)
where r(s), Q 0 (s), and Q τ (s) are polynomials with degrees satisfying the condition deg(r

) ≤ deg(Q 0 ) = deg(Q τ )
, and τ ≥ 0 is a constant. We first state the following result from [START_REF] Partington | H ∞ and BIBO stabilization of delay systems of neutral type[END_REF], which gives an explicit localization of the spectrum chain's asymptote for quasipolynomials of neutral systems.

Proposition 1 [START_REF] Partington | H ∞ and BIBO stabilization of delay systems of neutral type[END_REF] As such, under the condition |α| < 1, there are at most a finite number of poles of G(s) located in the right half plane. By scaling, this property extends to the half-plane to the right of any vertical line parallel to the imaginary axis, located to the right of the spectrum chain's asymptote. Note that for a retarded system, α = 0, and hence G(s) always has a finite number of poles in the right half plane. In contrast, for a neutral system, if |α| > 1, this property will vanish. Note further that in general, G(s) may contain multiple repeated poles, or equivalently, the quasipolynomial Q 0 (s) + Q τ (s)e -sτ may contain repeated zeros. In particular, in Boussaada and Niculescu (2016), it was shown that the algebraic multiplicity of a zero is bounded by the generic Polya/Szegö bound [START_REF] Pólya | Problems and Theorems in Analysis[END_REF], which is the degree of the quasipolynomial. Here by the degree of a quasipolynomial, we mean the sum of the number of the involved polynomials and their degrees, less one; for the quasipolynomial

Let α = lim |s|→∞ Q τ (s)/Q 0 (s).
Q 0 (s) + Q τ (s)e -sτ , its degree is deg(Q 0 ) + deg(Q τ ) + 1.
The Fritz John condition [START_REF] Bertsekas | Nonlinear Programming[END_REF]) is a firstorder necessary condition of optimality for constrained nonlinear programming problems whose objective function and constraints are differentiable, and the constraints may be equalities or inequalities. This class of problems can be stated as

min f (x) s.t. g i (x) ≤ 0, i ∈ I 1 = {1, 2, • • • , m} , h j (x) = 0, j ∈ I 2 = {m + 1, m + 2, • • • , k} , (2)
where f : R n → R, g i : R n → R and h j : R n → R all have continuous first-order partial derivatives in R n .

Proposition 2 The Fritz John Condition [START_REF] Bertsekas | Nonlinear Programming[END_REF]: If x * is an optimal solution of f (x) in (2), then there exists a row vector λ = [λ 0 , λ 1 , λ 2 , . . . , λ k ] such that:

λ 0 ∇f (x * ) + i∈I1 λ i ∇g i (x * ) + j∈I2 λ j ∇h j (x * ) = 0, (3) i∈I1 λ i g i (x * ) = 0, (4) λ i ≥ 0, i ∈ {0} ∪ I 1 , (5) λ = 0, ( 6
)
where ∇θ(x) = [∂θ(x)/∂x 1 , • • • , ∂θ(x)/∂x n ] denotes the gradient of θ(x).
It is useful to note that the Fritz John condition is applicable to problems where the constraints may not satisfy certain regularity conditions, specifically when the gradient vectors

∇g i (x * ) (i = 1, • • •, m) and ∇h * j (j = m + 1, • • •, k
) are linearly dependent, to which the seemingly better-known Karush-Kuhn-Tucker condition [START_REF] Luenberger | Linear and Nonlinear Programming[END_REF] fails to apply. On the other hand, when the gradient vectors are linearly independent, the Fritz John condition reduces to Karush-Kuhn-Tucker criterion. For a comprehensive treatment of Fritz John condition, we refer to [START_REF] Mangasarian | The Fritz John necessary optimality conditions in the presence of equality and inequality constraints[END_REF]; [START_REF] Bertsekas | Nonlinear Programming[END_REF].

Finally, we shall also need the concept of concavity, quasi-concavity, and pseudo-concavity.

Definition 3 [START_REF] Boyd | Convex Optimization[END_REF]) A function ϑ : R n → R is said to be concave over a convex set S ⊂ R n if for any x, y ∈ S,

ϑ((1-γ)x+γy) ≥ (1-γ)ϑ(x)+γϑ(y), γ ∈ [0, 1], (7)
and quasi-concave if

ϑ((1 -γ)x + γy) ≥ min{ϑ(x), ϑ(y)}, γ ∈ [0, 1]. (8)
Apparently, every concave function is quasi-concave.

Concave and quasi-concave functions admit a unique maximum, which can be found using convex optimization methods.

Definition 4 [START_REF] Boyd | Convex Optimization[END_REF]) A differentiable function ϑ : R n → R is said to be pseudo-concave over a convex set S ⊂ R n if for any x, y ∈ S,

∇ϑ(x) (y -x) ≤ 0 ⇒ ϑ(y) -ϑ(x) ≤ 0. ( 9
)
Pseudo-concave functions possess properties essentially similar to those of concave and quasi-concave functions.

Likewise, a unique maximum exists and can be computed in a similar manner. In particular, for a univariate pseudo-concave optimization problem, i.e., a unimodal pseudo-concave problem, the optimal solution can be solved efficiently by gradient-based and bisection approaches.

The delay margin and PID control

Consider the feedback system depicted in Figure 1, in which K(s) represents a finite-dimensional LTI controller, P τ (s) denotes the plant containing a constant but uncertain delay τ , with transfer function given by P τ (s) = P 0 (s)e -τ s , τ ≥ 0.

Here P 0 (s) represents a strictly proper finite dimensional delay-free plant. The controller K(s) of interest is one

   P s   K s r y Figure 1. A tracking control system.
of PID controllers, i.e., K(s) = K P ID (s),

K P ID (s) = k p + k i s + k d s. ( 10
)
The signals r and y represent the reference input and system output. Let

τ (k p , k i , k d ) = sup µ ≥ 0 : K P ID (s) stabilizes P τ (s), ∀τ ∈ [0, µ) .
Then τ (k p , k i , k d ) is the delay range within which the plant P τ (s) can be robustly stabilized by the PID controller K P ID (s) with the fixed gains (k p , k i , k d ). The delay margin of the system corresponds to the maximal range achievable by all PID controllers:

τP ID = sup τ (k p , k i , k d ) : K P ID (s) stabilizes P τ (s), ∀τ ∈ [0, µ) .
In this paper we consider the first-order delay plant

P τ (s) = 1 s -p e -τ s , τ ≥ 0, (11)
where p ≥ 0. For this class of plants, we have the following characterization of the delay margin, adapted from [START_REF] Ma | Delay margin of low-order systems achievable by PID controllers[END_REF].

Proposition 5 Let P τ (s) be given by ( 11). Then under the conditions k p > p, k i > 0, and

|k d | < 1, τ (k p , k i , k d ) = tan -1 ω0 p ω 0 + tan -1 k d ω0-k i ω 0 kp ω 0 , (12)
where ω 0 > 0 is given by

ω 2 0 = k 2 p -2k d ki-p 2 1-k 2 d + k 2 p -2k d ki-p 2 1-k 2 d 2 + 4 k 2 i 1-k 2 d 2 .
Furthermore,

τP ID = 2 p . (13)
Let the open-loop transfer function of the delay-free system be defined by

L 0 (s) = P 0 (s)K P ID (s). (14)
In the standard feedback design, the open-loop gain is to be shaped in such a way that |L 0 (jω)| 1 at low frequencies and |L 0 (jω)| 1 at high frequencies, thus necessitating the specification that |L 0 (0)| > 1 and |L 0 (∞)| < 1. Here the former loopshaping specification is required to ensure such performance objective as attenuation of low-frequency disturbances, and the latter is imposed for reduction of high-frequency noises and to maintain stability robustness against high-frequency modeled dynamics. With the integral control in presence, the condition |L 0 (0)| > 1 is rendered moot whenever k i > 0. The requirement |L 0 (∞)| < 1, on the other hand, coincides with the condition |k d | < 1. Note that for K P ID (s) to stabilize P 0 (s), additionally, it is necessary that k p > p, and k i > 0. In summary, we assume in the sequel that k p > p, k i > 0, and |k d | < 1.

Pole placement by PID control

While maintaining the delay robustness stipulated in Proposition 5, it is desirable for the PID controller to also achieve desirable transient behavior. In this section we provide a design based on multiplicity-induced dominancy, which enforces the closed-loop poles as dominant multiple poles, so that the closed-loop response has a prescribed decay rate; in other words, the rightmost poles of the closed-loop system are placed at a same location. For this purpose, we first note that with the first-order plant P τ (s) given by ( 11) and the PID controller given by ( 10), the closed-loop quasipolynomial is found as

∆(s) = Q 0 (s) + Q τ (s)e -sτ Q 0 (s) = s 2 -sp, Q τ (s) = k d s 2 + k p s + k i . ( 15
)
Our first main result of this paper is stated as follows.

Theorem 6 Let P τ (s) be given by ( 11)and K P ID (s) by ( 10). Then the following statements are true.

i) For arbitrary real parameters k p , k i , k d and arbitrary positive delay τ , the multiplicity of any root of the quasipolynomial ∆(s) is less than or equal to 5. ii) The quasipolynomial ∆(s) admits a multiple real root at

s ± = τ p -6 ± τ 2 p 2 + 12 2τ ( 16
)
with algebraic multiplicity greater than 3 if and only if

                   k d = (4 + 2 τ s ± -τ p) e τ s ± 2 , kp= - 8 τ + τ 2 s ± p-18-12 τ s ± e τ s ± τ , k i = (τ s ± + 3) τ 2 p 2 -(12τ s ± +60) τ p+108+84 τ s ± e τ s ± 2τ 2 . ( 17
)
iii) If ( 17) is satisfied, then under the condition τ < τP ID = 2/p, s = s + is the rightmost root of the quasipolynomial ∆(s) with multiplicity equal to 4.

PROOF. i)

The statement follows directly from the Pólya and Szegö bound (Pólya and Szegö, 1972), since the degree of the quasipolynomial ∆(s) is equal to 5.

ii) For a root s of the quasipolynomial ∆(s) given in ( 15), it follows that

e -τ s = -s 2 + ps k d s 2 + k p s + k i . ( 18
)
The root s has algebraic multiplicity greater than three if and only if in the ideal 18), the first-, second-, and third-order derivatives of ∆(s) are all zero at s. This results in three algebraic equations in 6 unknowns k i , k p , k d , τ, p, s:

I 4 =< ∂ s ∆, ∂ 2 s ∆, ∂ 3 s ∆ > generated by (
0 =τ k d s 4 + k p -k d p τ s 3 + k i -k p p τ + k d p + k p s 2 + (2 k i -τ pk i ) s -pk i , 0 = -τ 2 k d s 4 + -k p + k d p τ 2 + 4 τ k d s 3 + k p p -k i τ 2 + -4 k d p + 2 k p τ s 2 + τ 2 pk i + 2 k p + 2 k d p -2 τ k p p s + 2 k i , 0 =τ 3 k d s 4 + k p -k d p τ 3 -6 τ 2 k d s 3 + k i -k p p τ 3 + 6 k d p -3 k p τ 2 + 6 τ k d s 2 + -pk i τ 3 + 3 τ 2 k p p -6 τ k d p s (19)
By variable elimination, we obtain the following three sets of admissible solutions:

(1)

k i = k p = k d = 0. (2) k i = s = 0. (3) k p = -2 (τs±+8)τp-12 τ s±-18 k d τ (τ (2s±-p)+4) , k i = (τs±+3)τ 2 p 2 +(-12 τ s±-60)τ p+84 τ s±+108 k d (2 τ s±+4-τ p)τ 2 .
Evidently, the first two cases are not applicable and hence are discarded. In Case (3), substituting k p and k i into (18) yields the expression of k d in (17).

iii) The proof of the statement iii) is somewhat longwinding, and so for it to be more accessible, we divide the proof into three parts, each of which is preceded by a subheading.

iii-a) Establishing the multiplicity of s + To determine the exact multiplicity of s + , we consider the transformation s → (z/τ ) + s + and

∆(z) = τ 2 ∆ z τ + s + .
Denote ρ = τ 2 p 2 + 12. Under the condition 0 < τ < τP ID = 2/p, we have

2 √ 3 < ρ < 4. ( 20
)
A direct calculation gives rise to

∆(z) =z 2 + (ρ -6)z + (12 -3ρ) + ρ -2 2 z 2 + (2ρ -6)z + (3 ρ -12) e -z .
(21) It follows that s + is a root of ∆(s) with a certain algebraic multiplicity if and only if z = 0 is a root of ∆(z) with the same algebraic multiplicity. We then proceed to show that ∆(z) has a root at the origin with algebraic multiplicity 4. Toward this end, denote

Q0 (z) = z 2 + (ρ -6)z + (12 -3ρ) Qτ (z) = ρ -2 2 z 2 + (2ρ -6)z + (3 ρ -12)
and write

∆(z) = Q0 (z) + Qτ (z)e -z .
The quasipolynomial ∆(z) has a root on the imaginary axis if and only if ∆(jω) = 0. By a direct calculation, we find

| Q0 (jω)| 2 -| Qτ (jω)| 2 = 1 - ρ -2 2 2 ω 4 .
Since ρ < 4, it follows that for any ω = 0, | Q0 (jω)| > | Qτ (jω)|, and hence ∆(jω) = 0. To show that ∆(0) = 0, we obtain the Taylor series of ∆(z), given as

∆(z) = ρ 24 z 4 + o(|z| 4 ). ( 22
)
This enables us to conclude that ∆(z) has a root at the origin with multiplicity 4, and hence that ∆(s) has a root s + with multiplicity 4.

iii-b) Localizing the right half plane zeros of ∆(z)

We next show that ∆(z) has no root in the open right half plane, i.e., for all z such that Re(z) > 0. For this purpose, we first provide an auxiliary result, which is needed in the subsequent proof. This result is summarized in the following claim.

Claim Let z 0 = x 0 + j ω 0 , x 0 > 0, ω 0 ≥ 0 be a root of ∆(z). Then, ω 0 < π.

To establish this claim, we note that if z 0 is a root of ∆(z), i.e., ∆(z 0 ) = 0, then necessarily

| Q0 (x 0 + j ω 0 ) | 2 e 2 x0 =| Qτ (x 0 + j ω 0 ) | 2 . ( 23
)
Consider then the function

F (x, ω) =| Qτ (x + j ω) | 2 -(1 + 2 x) | Q0 (x + j ω) | 2 .
Since e 2 x > 1+2 x for any x > 0, we have F (x 0 , ω 0 ) > 0, for any right half plane root z 0 = x 0 + j ω 0 of ∆(z). The idea then is to show that if z 0 = x 0 + j ω 0 is a root of ∆(z), with x 0 > 0 and ω 0 ≥ π, then F (x 0 , ω 0 ) ≤ 0, thus resulting a contradiction. Toward this end, we set Ω = ω 2 and find the explicit expression of F (x, ω) as

F (x, √ Ω) = ρ 2 -4 ρ -8 x 4 Ω 2 + x 2 ρ 2 -12 ρ -8 x + 48 2 Ω -2 x 5 + (ρ -8) (ρ -12) 4 x 4
+ 24 (ρ -4)

x 3 + 18 (ρ -4) 2 x 2 .
Fix x > 0. The discriminant of the second-order polynomial in Ω is given by D(x) = x 2 D(x), where

D(x) = 16 ρ 2 -12 x 2 -24 (ρ -6) (ρ -4) 2 x -18 ρ (ρ -4) 3 .
Evidently, for ρ ∈ (2 √ 3, 4), D(x) > 0 for any x > 0. This means that the polynomial F (x,

√ Ω) has two real roots, which are given at the top of the next page.

Here Ω + (x) is the larger of the two roots. Accordingly, F (x, √ Ω) can be factorized as

F (x, √ Ω) = ρ 2 -4 ρ -8 x 4 Ω -Ω -(x) Ω -Ω + (x) .
(24) By judiciously setting ρ = 2 √ 3 or ρ = 4, we may establish the bound Ω + (x) ≤ Γ(x), where

Γ(x) = -x 2 -3 √ 3x + 15 x 2 + (-228 x + 468) √ 3 + 4 x 2 + 396 x -810.
Taking the derivative of Γ(x) gives rise to

Γ (x) = -2 x -3 √ 3 + 15 2 + -228 √ 3 + 8 x + 396 2 (-228 x + 468) √ 3 + 4 x 2 + 396 x -810
, which yields the unique positive stationary point x * ≈ 2.139487015. It thus follows that Ω + (x) ≤ Γ(x * ) ≈ 4.961051386 < π 2 . From (24), it is clear that for F (x, ω) > 0 for any x > 0, it is necessary that ω < π. As such, for F (x 0 , ω 0 ) > 0 with x 0 > 0, ω 0 ≥ 0, it is necessary that ω 0 < π, whereas the condition F (x 0 , ω 0 ) > 0 in turn is necessary for (23) to hold. We have thus established the claim.

Ω ± (x) = -8 x 2 + ρ 2 -12 ρ + 48 x ± 2 (16 ρ 2 -192) x 2 -24 (ρ -6) (ρ -4) 2 x -18 ρ (ρ -4) 3 x -ρ 2 + 4 ρ + 8 x .
iii-c) Establishing the dominance of s + To complete the proof, we employ the Fredholm integral representation of ∆(z) [START_REF] Kress | Linear Integral Equations[END_REF]. By a direct calculation, one can show that ∆(z) admits the Fredholm integral form

∆(z) = 1 0 K(t, z)q(t) dt, ( 25
)
where the kernel function K(t, z) and q(t) are given by

K(t, z) = z 4 e -tz 2 , q(t) = t(1 -t) (δ -4)t + 2 ,
respectively. Note that for ρ ∈ (2 √ 3, 4), q(t) is strictly positive over the interval (0, 1). Clearly, the Fredholm integral representation (25) can be simplified to

∆(z) = z 4 2 1 0 q(t) e -zt dt.
Assume that ∆(x 0 + jω 0 ) = 0 for some x 0 > 0 and ω 0 ≥ 0. Then necessarily,

         1 0 q(t)e -x0t cos(ω 0 t) dt = 0, 1 0 q(t)e -x0t sin(ω 0 t) dt = 0. ( 26
)
Since 0 < ω 0 < π, sin(ω 0 t) > 0 for t ∈ (0, 1). Likewise, q(t) > 0 for t ∈ (0, 1). Consequently,

1 0 q(t)e -x0t sin(ω 0 t) dt > 0,
which contradicts to (26). This proves that ∆(z) = 0 in the right half plane and hence establishes the dominance of s + , thus completing the proof.

Theorem 6 indicates that for any τ < τP ID = 2/p, by enforcing the PID control gains as in (17), we may always place the closed-loop dominant poles at + given in ( 16), which is a multiple negative real pole and hence guarantees that the closed-loop system is stable with a specified degree of stability. Note that for delay robustness (with respect to variation in delay) it is natural to impose the condition τ < τP ID = 2/p. Note also that implicit in the selection of the PID gains in accordance with (17), the derivative control gain k d satisfies the condition |k d | < 1 Under the condition τ < τP ID = 2/p, we have τ 2 p 2 + 12 < 4, and hence log |k d |/τ < s + . This allows us to conclude that ∆(s) has only a finite number of roots in C + s+ . Theorem 6 demonstrates further that ∆(s) in fact has no root in this region. Note additionally that the root s -is of no interest, since it leads to k d < -1; in other words, the resulting PID controller can not stabilize P 0 (s).

A number of remarks are now in order.

Remark 7 To further illustrate the proof of Theorem 6, we plot Ω + (x) and Γ(x) in Figure 2, for ρ ∈ (2 √ 3, 4) and a range of positive real x. One can see that within this range, Γ(x) serves as an upper bound of Ω + (x) independent of ρ, which does not exceed the value of 5. Figure 3 shows that Γ(x) has a unique maximum, and outside a certain range of x, it assumes negative values. Remark 8 The MID scheme can be employed for the analytic design and tuning of PI controllers, with dominating multiple closed-loop poles as well. Indeed, in the case of a PI controller, the closed-loop characteristic function is given by

∆(s) = s 2 -sp + k p s + k i e -sτ .
Define analogously the delay margin achievable by PI controllers, and denote it by τP I . Then from [START_REF] Ma | Delay margin of low-order systems achievable by PID controllers[END_REF], it is known that for the first-order plant (11), τP I = 1/p. By selecting s + = τ p -4 + τ 2 p 2 + 8 2τ

and

k p = τ 2 p 2 + 8 -2 e τ s+ τ , k i = (10 -pτ ) τ 2 p 2 + 8 + 2τ p -τ 2 p 2 -28 e τ s+
2τ 2 , one can show analogously that under the condition τ < τP I = 1/p, s = s + is the dominating root of ∆(s) with multiplicity 3.

We end our discussions by commenting briefly on the analytical MID tuning rule specified by (17). It is worth noting that the classical Ziegler-Nichols tuning rules, widely used for tuning PID controllers for stable processes, may have serious drawbacks and are not readily applicable to unstable plants [START_REF] Åström | Revisiting the Ziegler-Nichols step response method for PID control[END_REF]; for example, the Ziegler-Nichols' first and second methods rely on open-loop step and sinusoidal responses, which for open-loop unstable systems are difficult to measure. For this reason, PID controller tuning rules for unstable processes have been developed by focusing on the direct synthesis of the closed-loop response or the numerical optimization of a certain closed-loop performance objective, resulting in a wide variety of tuning rules that differ widely from each other (see, e.g., (O'Dwyer, 2013; Thomas and Ganesan, 2020) and the references therein). Unlike these earlier rules, the MID tuning rule presented herein adds an analytical tool to the available repertoire, which, as shown in the sequel, possesses desirable robustness properties.

Robustness of MID design

In this section we show that while tuned to the delay parameter and the value of the plant unstable pole, the PID controller ( 17) by the MID design can nevertheless render desirable robustness properties. In this vein, we first show that the MID design ensures a certain amount of robustness against uncertain delay, which in the limit approaches to its maximum, namely, the delay margin. This result is summarized in the following theorem.

Theorem 9 Let (k d , k p , k i ) be given in ( 17). Then the following statements are true.

i) For any τ < τP ID = 2/p,

     0 <k d < 1, k p > p, k i > 0. ( 27
)
The PID controller K P ID (s) stabilizes P τ (s) for all τ ∈ [0, τ ), where

τ = tan -1 ω0 p ω 0 + tan -1 k d ω0-k i ω 0 kp ω 0 ( 28
)
with ω 0 given by

ω 2 0 = k 2 p -2k d ki-p 2 1-k 2 d + k 2 p -2k d ki-p 2 1-k 2 d 2 + 4 k 2 i 1-k 2 d 2 . (29) ii) When τ → τP ID = 2/p,      k d → 1, k p → p, k i → 0, ( 30
)
and s + → 0. The PID controller K P ID (s) stabilizes P τ (s) for all τ ∈ [0, τP ID ). PROOF. To prove i), we first note that

τ s + = τ p -6 + (τ p) 2 + 12 2 . ( 31
)
For brevity, denote τ p = χ. Then for any τ < 2/p, χ < 2. Accordingly,

τ s + = χ -6 + χ 2 + 12 2 ,
which allows us to rewrite (17) as

           k d = f (χ)e τ s+ , k p = g(χ) τ e τ s+ , k i = h(χ) τ 2 e τ s+ .
where

                             f (χ) = χ 2 + 12 2 -1, g(χ) = 6 - χ 2 χ 2 + 12 - χ 2 2 + χ -18, h(χ) = χ 2 2 -6 χ + 42 χ 2 + 12 + χ 3 2 -6 χ 2 + 18 χ -144.
Obviously, f (χ) > 0. For any χ < 2, we have f (χ) < 1. Since e τ s+ < 1, k d = f (χ)e τ s+ < 1. Consider next the function ĝ(χ) = g(χ)e g1(χ) -χ,

where

g 1 (χ) = χ -6 + χ 2 + 12 2 .
We prove below that ĝ(χ) is monotonically increasing on [0, 2). For this purpose we consider its first-and secondorder derivatives χ) .

ĝ (χ) = (χ -5) χ 2 + 12 -χ 2 + 5 χ -8 2 e g1(χ) -1 and ĝ (χ) = -χ 2 + 3 χ -3 2 + χ 3 -3 χ 2 + 9 χ -18 2 χ 2 + 12 e g1(
One can easily show that ĝ (χ) is negative since both the polynomials g 2 (χ) = -χ 2 + 3 χ -3 and g 3 (χ) = χ 3 -3 χ 2 + 9 χ -18 are negative on the interval [0, 2). In other words, ĝ (χ) is a decreasing function over the interval [0, 2), which in turn means that inf

χ∈[0,2) ĝ (χ) = ĝ (2) = 0.
As such, on [0, 2), ĝ (χ) > 0, which implies that ĝ(χ) is monotonically increasing, and for 0 < χ < 2, ĝ(χ) = g(χ)e g1(χ) -χ > ĝ(0) > 0, or equivalently g(χ)e g1(χ) > χ. Substituting now χ = τ p, we have

k p = g(χ) τ e τ s+ > p.
Finally, consider the parameter k i . We show that h(χ) > 0. Toward this end, we first find

h (χ) = 3 (χ -2) (χ -6) χ 2 + 12 + χ + 24 2 χ 2 + 12 . For χ < 2, it is immediate that (χ -6) χ 2 + 12 + χ + 24 < 0,
and hence h (χ) < 0. In other words, h(χ) is monotonically decreasing on [0, 2), and inf χ∈[0,2) ĥ(χ) = ĥ(2) = 0.

This proves that k i > 0. Under the condition ( 27), it follows from Proposition 5 that P τ (s) can be stabilized robustly by K P ID (s) for all τ ∈ [0, τ ), where τ is given by ( 28). This completes the proof for the statement i).

To prove the statement ii), it suffices to note that when τ → τP ID , s + → 0.

Remark 10 From ( 16), a direct calculation shows that

ds + dτ = 3 τ 2 1 - 2 τ 2 p 2 + 12 .
Hence, for any τ < 2/p, ds + /dτ > 0, that is, s + increases monotonically with τ , thus demonstrating a tradeoff between the delay τ selected and the relative stability of the closed-loop system. It is intuitively plausible to contend that an inherent tradeoff may exist as well between the dominant root s + and the delay margin τ . The relationship between these two quantities, however, appears highly complex to prevent from an analytical characterization, whilst examples can be constructed to show that a monotone relationship does exist (cf. Section 6). Interestingly, the example presented in Section 6 also shows that not only does the MID design guarantee a certain amount of delay robustness, but it in fact yields a delay margin, surprisingly, greater than the value of τ selected.

It is of interest to analyze the robustness of the MID design to the variation in the plant unstable pole. For this purpose, we examine the uncertain plant

Pτ (s) = 1 s - p e -τ s , τ ≥ 0, (32)
with an uncertain pole p ∈ [p -ν, p + ν], for some prescribed ν > 0. Define analogously

∆(s, p) = s 2 -sp + (k d s 2 + k p s + k i )e -sτ .
Under the condition τ < τP ID = 2/p, let (k p , k i , k d ) be given by the MID design ( 17) with the nominal value p of the plant pole. We first analyze the variation of the dominant pole s + under the perturbation p -p, which can be formulated as an eigenvalue perturbation problem. Specifically, by exploiting the results of [START_REF] Chen | Stability analysis of polynomially dependent systems by eigenvalue perturbation[END_REF], we assert that the multiple pole s + splits into four branches, which in the neighborhood of s + can be described asymptotically as

s (k) + = s + +   -4! ∂∆(s+, p) ∂ p | p=p ∂ 4 ∆(s, p) ∂s 4 | s=s+   1/4 (p -p) 1/4 +o((p -p) 1/4 ), k = 1, 2, 3, 4.
Here for a complex number z, we denote its m-th roots by

z 1/m = |z| 1/m e j(2(k-1)π+ z)/m , k = 1, 2, • • • , m.
It is evident that

∂∆(s + , p) ∂ p | p=p = -s + .
Furthermore, since ∆(s, p) = ∆(s) = ∆(τ (s -s + ))/τ 2 , it follows from ( 22) that

∂ 4 ∆(s, p) ∂s 4 | s=s+ = ρτ 2 .
Consequently, we find that for k = 1, 2, 3, 4,

s (k) + = s + + 4! ρτ 2 1/4 (s + ) 1/4 (p -p) 1/4 + o((p -p) 1/4 ).
(33) Note that for p > p, since s + < 0, the dominant pole s + will split into four branches with asymptotic departure angles as π/4, 3π/4, 5π/4, and 7π/4, respectively. As p increases, two of the four poles at s + will emerge as a pair of dominate conjugate complex poles with asymptotes oriented at the angles of π/4 and -π/4.

Tradeoff between delay robustness and steady-state tracking

With its integral control action, PID control generally performs well with set-point reference tracking. On the other hand, Theorem 9 suggests that to gain a maximal delay margin, the integral control effort must be reduced, thus demonstrating a tradeoff between steadystate tracking and delay robustness, which more generally can be attributed as a manifestation of the tradeoff between control performance and system robustness.

In this section, we characterize this tradeoff in a precise manner.

For the feedback configuration shown in Figure 1, the tracking error signal is defined by

e(t) = y(t) -r(t).
Let the Laplace transforms of the reference input r(t) and the tracking error signal e(t) be denoted as r(s) and ê(s). It follows that

ê(s) = 1 1 + L 0 (s)e -τ s r(s), (34)
where L 0 (s) is given by ( 14). Provided that the closedloop system is stable, the steady-state tracking error is determined as

e(∞) = lim t→∞ e(t) = lim s→0 sê(s).
Evidently, with a step reference signal, perfect steadystate tracking is achieved whenever the integral control is applied, that is, the steady-state tracking error is zero.

Consider then the more demanding task of tracking a ramp reference signal:

r(t) =    t t ≥ 0 0 t < 0.
In this case, the integral control is clearly necessary; otherwise, the asymptotic tracking error would be infinitely large. For the first-order plant given in (11), the absolute tracking error becomes

|e(∞)| = lim s→0 s 1 + L 0 (s)e -τ s r(s) = p k i ,
whenever the closed-loop system is stable. In other words, to ensure a small |e(∞)|, it is necessary to make k i large. To meet this requirement, we enforce the specification given in

Assumption 11 k i ≥ σ > 0.
Under this assumption, define

Ξ(σ) = (k p , k i , k d ) : k p > p, |k d | < 1, k i ≥ σ , and τP ID (σ) = sup τ (k p , k i , k d ) : (k p , k i , k d ) ∈ Ξ(σ) .
Then τP ID (σ) defines the maximal delay margin achievable contingent upon the tracking performance desired. Note that τP ID (0) = τP ID . The following result characterizes this delay margin, and henceforth the tradeoff in achieving delay robustness and the tracking performance.

Theorem 12 Let P τ (s) be given by (11). Then the following statements hold:

i) τP ID (σ) = sup p<kp< √ p 2 +2σ τ (k p ), (35)
where τ (k p ) = lim

k d →1 τ (k p , σ, k d ), (36)
with τ (k p , k i , k d ) = τ in (28), and ω 0 given by

ω 0 = σ p 2 -k 2 p + 2σ . ii) τ (k p ) is pseudo-concave for p < k p < p 2 + 2σ.
iii) τP ID (σ) is strictly decreasing with σ ≥ 0.

iv)

τP ID (σ) ≤    tan -1 √ σ 2 p √ σ 2 p    2 p . (37)
It is clear from i) that for a first-order unstable system, the optimal derivative control coefficient k d lies (asymptotically) on the boundary of its allowable range, i.e., k d = 1, pointing to the fact that the derivative control has a positive effect to compensate for the time delay. On the contrary, the delay margin τP ID (σ) is achieved when the integral control gain is minimal, i.e., k i = σ, indicating that the integral control contributes no effect to enlarge the delay margin. The pseudo-concavity of τ (k p ) reveals that the exact determination of the delay margin amounts to solving a pseudo-concave unimodal problem, which can be solved using convex optimization or gradient-based methods. The monotonicity property established iii) demonstrates the tradeoff between delay robustness and tracking accuracy, which can be further observed by the explicit a priori bound stated iv).

We now provide the proof of Theorem 12. Similarly, we also provide subheadings as a guide to the proof.

PROOF. i-a) Reformulation as a nonlinear programming problem We begin by recognizing the fact that the set Ξ(σ) parameterizes the PID coefficients such that the closed-loop system is stable, and the open-loop gain satisfies the roll-off condition |L 0 (∞)| < 1, together with the steady-state tracking specification prescribed by Assumption 11. Under these constraints, we seek to maximize τ (k p , k i , k d ). Our first step is to recast the maximization problem as one of constrained nonlinear programming:

min f (k p , k i , k d , ω 0 ) s.t. g 1 (k p , k i , k d , ω 0 ) = p -k p ≤ 0, g 2 (k p , k i , k d , ω 0 ) = k d -1 ≤ 0, g 3 (k p , k i , k d , ω 0 ) = -k d -1 ≤ 0, g 4 (k p , k i , k d , ω 0 ) = σ -k i ≤ 0, h(k p , k i , k d , ω 0 ) = (1 -k 2 d )ω 4 0 + (p 2 -k 2 p + 2k d k i )ω 2 0 -k 2 i = 0.
where

f (k p , k i , k d , ω 0 ) = - tan -1 ω0 p ω 0 - tan -1 k d ω0-k i ω 0 kp ω 0 .
Note that in this formulation, the first four inequality constraints represent the set Ξ(σ), where the strict inequalities are modified to enable the equivalence between the infimum and minimum. The last equality constraint characterizes the crossover frequency ω 0 , i.e., the solution to h(k p , k i , k d , ω 0 ) = 0. It is clear that

τP ID (σ) = sup (kp,ki,k d )∈Ξ(σ) τ (k p , k i , k d ) = min g1,g2,g3,g4,h f (k p , k i , k d , ω 0 ). i-b) Establishing the necessity of k i = σ, k d = 1
for optimality Note that since the equation h(k p , k i , k d , ω 0 ) = 0 admits a unique solution given by ( 29), the objective function f (k p , k i , k d , ω 0 ) is differentiable. This allows us to invoke the Fritz John condition, which gives rise the first-order conditions

λ 0 f (k * p , k * i , k * d , ω * 0 ) + 4 i=1 λ i g i (k * p , k * i , k * d , ω * 0 ) + λ 5 h(k * p , k * i , k * d , ω * 0 ) = 0, 4 i=1 λ i g i (k * p , k * i , k * d , ω * 0 ) = 0, λ i ≥ 0, i ∈ {0, . . . , 4}, λ = [λ 0 , λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ] = 0.
We then examine the partial derivatives of f (k p , k i , k d , ω 0 ), h(k p , k i , k d , ω 0 ) and g i (k p , k i , k d , ω 0 ), i = 1, . . . , 4 with respect to k d , k i and k p , which results in the equations

∂ ∂k d : λ 0 1 1 + Π 2 1 k * p -λ 2 +λ 3 +2λ 5 (k * d ω * 2 0 -k * i )ω * 2 0 = 0, (38) ∂ ∂k i : λ 0 1 1 + Π 2 1 k * p -λ 4 ω * 2 0 + 2λ 5 (k * d ω * 2 0 -k * i )ω * 2 0 = 0, (39)
∂ ∂k p : λ 0 1 1 + Π 2 k * d - k * i ω * 2 0 1 k * 2 p -λ 1 -2λ 5 k * p ω * 2 0 = 0, (40)
where

Π = (k * d ω * 0 -k * i /ω * 0 )/k * p .
In light of the condition

4 i=1 λ i g i (k * p , k * i , k * d , ω * 0 ) = 0, it follows that g i (k * p , k * i , k * d , ω * 0 )
= 0 whenever λ i > 0, which means that the constraint g i ≤ 0 is active. Note that the equations ( 38) and ( 39) give rise to

-λ 2 + λ 3 = -λ 4 ω * 2 0 . (41) Since g 2 (k * p , k * i , k * d , ω * 0 ) = 0 and g 3 (k * p , k * i , k * d , ω * 0 )
= 0 are mutually exclusive, λ 2 and λ 3 cannot be both positive. Hence for (41) to hold, either λ 2 = λ 3 = λ 4 = 0 or λ 3 = 0, λ 2 = λ 4 ω * 2 0 . Consider first the case λ 2 = λ 3 = λ 4 = 0. The discussion of this case can be further divided into two parts.

1) k *

d ω * 2 0 -k * i < 0. If λ 0 = 0, then λ 5 = 0 by (38). Since λ = 0, we have λ 1 = 0. This, however, contradicts to (40). Hence, we assert that λ 0 > 0. Under this condition, it follows from (38) that λ 5 > 0. Nevertheless, that both λ 0 > 0 and λ 5 > 0 contradicts to (40).

2) k * d ω * 2 0 -k * i ≥ 0. If λ 0 = 0, then from (40) we have λ 1 + 2λ 5 k * p ω * 2 0 = 0, which in turn means that λ 1 > 0 and λ 5 < 0; otherwise, it is necessary that λ = 0, which is a contradiction. Similarly, if λ 0 > 0, then from (38) we have λ 5 < 0. This again allows us to conclude from (40) that λ 1 > 0. Hence in either case, since λ 1 > 0, the constraint g

1 (k p , k i , k d , ω 0 ) ≤ 0 is active, i.e., g 1 (k * p , k * i , k * d , ω * 0 ) = p -k * p = 0.
Likewise, that λ 5 < 0 implies that the equality constraint h(k p , k i , k d , ω 0 ) = 0 is active. Solve these two equations yields

(1 -k * d )ω * 2 0 = -k * i ,
which clearly is a contradiction. In summary, we have shown that the solution λ 2 = λ 3 = λ 4 = 0 invariably leads to a contradiction and hence is not possible.

Consider then the case λ 3 = 0, λ 2 = λ 4 ω * 2 0 . Clearly, neither λ 2 nor λ 4 can be equal to zero. In other words, λ 2 > 0, and λ 4 > 0. As a result, the constraints g 2 (k p , k i , k d , ω 0 ) ≤ 0 and g 4 (k p , k i , k d , ω 0 ) ≤ 0 are active and in turn k

* d = 1, k * i = σ. With k * d = 1 and k * i = σ, we obtain from the constraint h(k * p , k * i , k * d , ω * 0 ) = 0 the solution ω * 2 0 = σ 2 p 2 -k * 2 p + 2σ
.

This establishes (36). Note that for k p ∈ (p, p 2 + 2σ), ω 0 ∈ ( σ/2, ∞).

ii) Proving the pseudo-concavity of τ (k p ) Denote by τ (k p ) = f (g(k p )), where

ω 0 = g(k p ) = σ p 2 -k 2 p + 2σ and f (ω 0 ) = tan -1 ω0 p ω 0 + tan -1 ω 2 0 -σ √ (p 2 +2σ)ω 2 0 -σ 2 ω 0 . ( 42
) It follows that τ (k p ) = f (ω 0 )g (k p ), with f (ω 0 ) = ω0 p 1 + ( ω0 p ) 2 -tan -1 ω 0 p -tan -1 φ(ω 0 ) + φ (ω 0 )ω 0 1 + φ 2 (ω 0 ) 1 ω 2 0 , where
φ(ω 0 ) = ω 2 0 -σ (p 2 + 2σ)ω 2 0 -σ 2 and φ (ω 0 ) = (p 2 + 2σ)ω 3 0 + p 2 σω 0 (p 2 + 2σ)ω 2 0 -σ 2 3 2 .
Clearly, g (k p ) > 0. We prove below that a unique kp exists such that τ (k p ) > 0 for p < k p < kp and τ (k p ) < 0 for kp < k p < p 2 + 2σ, so as to prove the pseudo-concavity of τ (k p ). Consider the function δ(ω 0 ) = ω 2 0 f (ω 0 ). We have δ(ω 0 ) → -π < 0 when ω 0 → ∞. On the other hand, when ω 0 → σ/2, we have tan -1 ω 0 p → -tan -1 φ(ω 0 ), which suggests that δ(ω 0 ) > 0 as ω 0 → σ/2. By the continuity of δ(ω 0 ), there exists at least one ω0 > σ/2 such that δ(ω 0 ) = 0; correspondingly, there exists some p < kp < p 2 + 2σ such that τ ( kp ) = 0. Next, we prove that this kp is unique, by showing that δ(ω 0 ) is a monotonically decreasing function. Indeed, taking the derivative of δ(ω 0 ) yields

δ (ω 0 ) = - 2 p ( ω0 p ) 2 (1 + ( ω0 p ) 2 ) 2 - ω 0 θ(ω 0 ) 1 + φ 2 (ω 0 ) 2 ,
where

θ(ω 0 ) = 2φ 2 (ω 0 )φ(ω 0 ) -(1 + φ 2 (ω 0 ))φ (ω 0 ), and
φ (ω 0 ) = - (2p 2 + 3σ)(p 2 + 2σ)σω 2 0 + p 2 σ 3 (p 2 + 2σ)ω 2 0 -σ 2 5 2 . Rewrite θ(ω 0 ) = ω 4 0 γ(ω 0 ) (p 2 + 2σ)ω 2 0 -σ 2 7 2
, where

γ(ω 0 ) = 2(p 2 + 2σ) 2 ω 4 0 + (p 2 + 2σ)(4p 2 σ + 7σ 2 )ω 2 0 + (6p 4 σ 2 + 3p 2 σ 3 + 2p 6 σ).
Obviously, γ(ω 0 ) > 0, and so θ(ω 0 ) > 0. As such, δ (ω 0 ) < 0 and δ(ω 0 ) is monotonically decreasing with ω 0 , and that the solution ω0 to the equation δ(ω 0 ) = 0 is unique. Since g (k p ) > 0 and τ (k p ) = δ(ω 0 )g (k p )/ω 2 0 , the uniqueness of kp is ascertained. We thus concludes that τ (k p ) is pseudo-concave.

iii) Proving the monotonicity of τP ID (σ) Let σ 1 < σ 2 . From the above proof, τP ID (σ 1 ) is uniquely achieved on the boundary k

* i = σ 1 for k i in the interval [σ 1 , ∞). Since for any σ 2 ∈ (σ 1 , ∞), i.e., σ 2 lies in the interior of [σ 1 , ∞), sup τ (k p , k i , k d ) : (k p , k i , k d ) ∈ Ξ(σ 1 ) > sup τ (k p , k i , k d ) : (k p , k i , k d ) ∈ Ξ(σ 2 ) . Thus, τP ID (σ 1 ) > τP ID (σ 2 ), that is, τP ID (σ) is strictly decreasing with σ > 0.
iv) The upper bound on τP ID (σ) In view of (42), since for ω 2 0 ≥ σ,

tan -1 ω 0 -σ ω0 k p ≤ tan -1 ω 0 k p ≤ tan -1 ω 0 p , it follows that τ (k p ) ≤ 2 tan -1 ω0 p ω 0 .
Furthermore, since tan -1 ω 0 p ω0 is monotonically decreasing with ω 0 , we have sup

ω 2 0 ≥σ τ (k p ) ≤   tan -1 √ σ p √ σ p   2 p . For σ/2 < ω 2 0 < σ, since tan -1 (ω 0 -σ/ω 0 )/k p < 0, we have τ (k p ) ≤ tan -1 ω0 p ω 0 .
The maximum of this upper bound is achieved at ω 0 = σ/2, which leads to

sup σ/2<ω 2 0 <σ τ (k p ) ≤    tan -1 √ σ 2 p √ σ 2 p    2 p .
Taking the larger of the two upper bounds completes the proof.

An illustrative example

To illustrate the preceding results, we consider the delay plant (11) with p = 1. For this plant, the maximal delay margin is τP ID = 2. Select τ = 1. Tuning the PID controller gains as prescribed in (17) gives

k d = -1 + √ 13 2 e -5 2 + √ 13 2 , k i = - 263 4 + 73 √ 13 4 e -5 2 + √ 13 2 , k p = - 35 2 - √ 13 2 e -5 2 + √ 13 2 , (43)
which results in the quadruple real closed-loop pole s + = (-5 + √ 13)/2. Obviously, s + is a stable pole. Figure 4 exhibits the closed-loop spectrum distribution with the PID control gains given by (43). Figure 5 illustrates the relationship between τ selected in the PID design and the dominant pole s + , which serves to show the tradeoff between the delay tolerable and the transient performance achieved, such as the decay rate of the closed-loop system. Note that different values of p > 0 are examined. It is worth noting that the two quantities possess a monotone relationship. Figure 6 is helpful for assessing delay robustness, where different values of p > 0 are also selected. Interestingly, in all the cases, and for any τ selected (equivalently, s + is selected), the delay robustness range is greater than τ , which means that with the MID design, one can guarantee at least the same amount of robustness range as τ . Likewise, as τ increases, the dominant pole s + becomes closer to the imaginary axis, but the delay robustness range τ is increased. We next analyze the delay margin vs. steady-state tracking performance. Select p = 1.5 and σ = 1. Figure 7 demonstrates that τ (k p ) is indeed a pseudo-concave function. Figure 8 plots the delay margin σP ID (σ) for σ = 0, 0.2, 0.5 and 1, with p ∈ [0, 3]. It is clear that as σ increases, the delay margin decreases. For different values of σ, the reduction of the delay margin, due to steady-state tracking specification, can be significant.

Conclusion

Pole placement of time-delay systems is a challenging problem which often defies accurate analysis and design. The difficulty becomes exacerbated when using PID control, due to its restricted degree of design freedom. Additionally, complexity may also arise for delay plants with relative degree less than two, since PID control will result in a neutral closed-loop delay system which may have infinitely many poles in the right half plane. This paper presents a case study where the difficulties are judiciously mitigated. We developed a MID design and tuning scheme for first-order unstable delay plants, which has been shown to possess desirable properties. Specifically, we show analytically that with this scheme, a real dominant closed-loop pole can be assigned precisely, thus guaranteeing the exponential decay rate of the system response. While achieving such transient performance, the MID design is also shown to be capable of maintaining a certain level of stability robustness against an uncertain delay. The tradeoff of the design, versus that of achieving delay margin and steady-state tracking performance, is analyzed. The latter is shown to be computable by solving a unimodal convex optimization problem.

Extensions of this work may be pursued in several aspects. First, it is of immediate interest to contemplate on the problem of assigning dominant complex conjugate poles. This will further enrich the transient behavior of the closed-loop system, in gaining, e.g., a faster rise time and a shorter setting time. Likewise, it is of significant interest to seek extensions to second-and higherorder systems. Moreover, in practice, the implementation of the derivative control is often accompanied by a low-pass filter, resulting in the so-called filtered derivative control [START_REF] Åström | PID Controllers: Theory, Design, and Tuning[END_REF]. It is of interest to develop a companion MID design and tuning rule for such practical implementation. For these purposes, it appears that the techniques employed herein, which draw heavily upon the theories of algebraic geometry and nonlinear programming, need to be generalized. Finally, it will be of interest to explore applications of the MID method to practical design problems. In this vein, we note that the MID method is in general applicable to linear time-invariant systems in which time delays may exist in the plant dynamics or control action is implemented via delayed feedback. Successful applications of this method have been reported in, e.g., the control of drilling systems [START_REF] Auriol | Torsional-vibrations Damping in Drilling Systems: Multiplicity-Induced-Dominancy-Based Design[END_REF].

i)

  If |α| < 1 then the poles of G(s) of large modulii are asymptotic to the vertical line Re(s) ≈ log(|α|)/τ in the left half plane. The number of poles of G(s) to the right of Re(s) = log(|α|)/τ + is finite for any > 0. ii) If |α| > 1 then G(s) has infinitely many unstable poles, asymptotic to the vertical line Re(s) ≈ log(|α|)/τ in the right half plane.

Figure 2 .

 2 Figure 2. 3D-plot of Ω+(x) (yellow) and Γ(x) (blue). (cf. Theorem 9). As noted at the end of Section 2, this condition coincides with the open-loop roll-off requirement. Inadvertently, it also guarantees that only a finite number of roots of ∆(s) exist in the shifted right-half complex plane C + s+ = {s ∈ C, s.t. Re(s) > s + }. To see this point, we note that α = lim |s|→∞ Q τ (s)/Q 0 (s) = k d .

Figure 3 .

 3 Figure 3. The function Γ(x).

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. The distribution of the closed-loop spectrum with (kp, ki, k d ) given by (43), corresponding to p = τ = 1. The dominancy of s+ = -5- √ 13 2 is emphasized.

Figure 7 .

 7 Figure 7. Pseudo-concavity of τ (kp).

Figure 8 .

 8 Figure 8. Delay margin vs. tracking error.

Acknowledgments

The authors are partially financially supported by NSFC under Grant 61973060, the Hong Kong RGC under Project CityU 11203120, a public grant overseen by the French National Research Agency (ANR) as part of the "Investissement d'Avenir" program, through the "iCODE Institute project" funded by the IDEX Paris-Saclay.


Marco Carozzi
email: marco.carozzi@unimi.it

Mattia Rossana Monica Ferrara

Mattia Fumagalli

Mattia Sanna

Marcello Chiodini

Ales-Sia Perego

Alessandro Chierichetti

Stefano Brenna

Gianfranco Rana

Marco Acutis

Ales- Sia Perego

Field-scale ammonia emissions from surface spreading of dairy slurry in Po Valley

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

eutrophication of water with loss of biodiversity, respiratory diseases and the long-range transport of sulphur (S) and N [START_REF] Sutton | The exchange of atmospheric ammonia with vegetated surfaces. II: fertilised vegetation[END_REF][START_REF] Asman | Ammonia: emission, atmospheric transport and deposition[END_REF][START_REF] Erisman | Instrument development and application in studies and monitoring of ambient ammonia[END_REF][START_REF] Harper | Ammonia: measurement issues[END_REF]. Moreover, by 2020 NH 3 is expected to be the largest single contributor to acidification, eutrophication and formation of secondary particulate matter [START_REF] Ammann | Baseline Scenarios for the Clean Air for Europe (CAFE) Programme[END_REF]. The need of reliable NH 3 measurements at fieldscale becomes decisive (i) to promote abatement strategies, (ii) to derive emission factors to be used in national and international emission inventories, (iii) to validate models, (iv) to evaluate the ammonia exchange over natural surfaces in the continuum soil-plants-atmosphere domain, (v) to quantify the value of agronomic N-fertilisers. The ammoniacal losses from agriculture contributes to over 90% in Europe (EEA, 2011), where Po Valley (Northern Italy) is one of the most emitting region of the whole area [START_REF] Skjøth | Spatial and temporal variations in ammonia emissions -a freely accessible model code for Europe[END_REF]. Nevertheless, the lack of measured data at

INTRODUCTION

Agriculture is the primary source of gaseous ammonia (NH 3 ) in atmosphere and its emissions are mainly originated from the field application of animal manure and fertilisers (Gènermont et al., 1998;Sommer et al., 2001;[START_REF] Asman | Emissions of ammonia[END_REF], animal waste and grazing [START_REF] Jarvis | Ammonia volatilisation from agricultural land. The fertiliser society proceedings[END_REF]. Ambient NH 3 assumes an important role and growing interest among different atmospheric nitrogen reactive species as a key to mitigate the impact of nitrogen (N) on terrestrial ecosystems [START_REF] Sutton | Scope and Overview of the UNECE Expert Workshop on Ammonia[END_REF]. The environmental issues due to NH 3 emissions include mainly acidification of soils, Abstract: Po Valley (Northern Italy) is one of the major ammonia (NH 3 ) emitting regions of Europe, where the slurry spreading causes high NH 3 volatilisation, reducing its agronomic value and becoming a potential cause of environmental concerns. In autumn 2011 a field trial was conducted to estimate the NH 3 losses from the application of dairy slurry at rate of 57 m 3 ha -1 on bare soil. The emissions were estimated from surface application of dairy slurry by using an inverse dispersion modelling technique associated with long term exposure passive samplers and the measure of the atmospheric turbulence. NH 3 emissions levels resulted high within the first 24 hours from the spreading, reaching the 73% of the entire losses, with a maximum value of 163 μg m -2 s -1 after 3 hours and 20 minutes, whereas the 50% of the emissions was achieved after 10 hours. The phenomenon stopped after 168 hours with a total NH 3 losses equal to 44% of the total ammoniacal nitrogen (TAN) applied. Results showed and confirmed that surface application involves high NH 3 emissions and then alternative low-emission techniques have to be adopted. Keywords: ammonia emissions, inverse dispersion modeling, surface slurry application, passive samplers.

Riassunto:

In pianura Padana la distribuzione superficiale dei reflui zootecnici determina elevate emissioni di ammoniaca (NH 3 ), riducendo il loro valore agronomico e causando problemi ambientali. L'obiettivo di questo studio è la stima dell'emissione di NH 3 a seguito di una distribuzione superficiale di 57 m 3 ha -1 di reflui zootecnici su suolo nudo nell'autunno 2011. La stima è stata effettuata mediante l'applicazione di un modello per la dispersione degli inquinanti, associato alla misura della concentrazione dell'NH 3 mediante esposizione in pieno campo di campionatori a diffusione passiva e all'utilizzo di un anemometro sonico per la misura della turbolenza atmosferica. L'emissione di ammoniaca è stata elevata nelle prime 24 h dalla distribuzione, evidenziando un picco massimo dopo 3 ore e 20 minuti (163 μg m -2 s -1 ), laddove il 50% dell'intera emissione è stato raggiunto già a 10 ore dalla distribuzione. L'emissione di NH 3 è stata pari al 44% del totale di azoto ammoniacale applicato. I risultati mostrano e confermano come la distribuzione superficiale sia un metodo che determina alte perdite di ammoniaca e che quindi deve essere incentivato l'uso di tecniche alternative. Parole chiave: emissioni di ammoniaca, modelli a dispersione, distribuzione superficiale dei reflui, campionatori passivi.

characterized by silty-clay soil (Hypercalcic Calcisol). Maize crop was previously harvested in September and stubbles were chopped and left on the surface. A uniform dairy slurry application rate of 57 m 3 ha -1 was applied on 10 th of October using a 20 m 3 tank equipped with a splash plate. The spreading started from the longitudinal row passing to the centre of the field, close to measuring devices, towards the upwind edge, and then from the centre to the downwind side of the field. The slurry application (started at 8.15 a.m. and lasted 4 hours) supplied 68 kg N ha -1 of total ammoniacal nitrogen (TAN = NH 4 + + NH 3 ). The TAN content was the 63% of the total slurry N content, while the values of the dry matter and the pH were 30 g kg -1 and 7.5, respectively. Meteorological data were collected by a standard weather station, located close to the field. During the trial air temperature ranged from 0 to 25°C (mean value: 12.3°C), relative humidity was 29 to 98% (mean value: 65.5%). The mean of wind speed was 1.2 m s -1 with a maximum value of 4.7 m s -1 , whereas the main wind direction was SW. No rain events occurred in the sampling period. During the experiment ammoniacal (NH 4 -N) and nitrate (NO 3 -N) nitrogen, pH and water content were daily measured at 0-5 cm and 5-15 cm soil depth. Particularly, concentrations of soil NH 4 -N and NO 3 -N were performed with a KCl extraction and determined by spectrometric detection (FIAstar 5000 Analyzer, Foss Tecator, Denmark). These analysis were in agreement with the ISO 11732 (1997) andISO 13395 (1996) procedures. The values of soil pH were obtained in water solution with a soil-to-solution ratio of 1:2.5 (weight/volume), whereas the soil water content (SWC) was determined gravimetrically for each soil sample by the oven-drying method.

Air ammonia measurements

The air NH 3 concentration was quantified through the exposure of the passive samplers ALPHA (Adapted Low-cost Passive High Absorption) developed by [START_REF] Tang | Development and Types of Passive Samplers for Monitoring Atmospheric NO 2 and NH 3 Concentrations[END_REF] and Sutton et al. (2001a). The operating principle of ALPHA samplers is the capture of gaseous NH 3 on acid support coated with citric acid. These tools are designed to measure NH 3 air concentration less than 1 μg m -3 [START_REF] Leith | Quantifying dry NH 3 deposition to an ombrotrophic bog from an automated NH 3 field release system[END_REF] to over 4 mg m -3 [START_REF] Carozzi | A note on estimating the effect of a limited fetch on micrometeorological evaporation measurements[END_REF]. Samplers were placed both in the centre of the field to measure the NH 3 concentration from the slurry application (C) and 1 km away from the field and from any other known source of NH 3.

field scale in such region [START_REF] Valli | national inventory of ammonia and greenhouse gas emissions from agriculture in Italy[END_REF], was only recently filled [START_REF] Carozzi | A note on estimating the effect of a limited fetch on micrometeorological evaporation measurements[END_REF]. NH 3 losses from field-applied manure, particularly slurry, were measured in many European experiments [START_REF] Søgaard | Ammonia volatilization from fieldapplied animal slurry -the ALFAM model[END_REF]Sintermann et al., 2011a). Due to the sticky characteristics given by its polar configuration, NH 3 is capable to bind and to be released from solid surfaces, resulting in biased measurement of the emission. However, despite many techniques have been developed [START_REF] Brodeur | Technical note: laboratory evaluation of a tunable diode laser system for eddy covariance measurements of ammonia flux[END_REF], a standardized method is not yet available. The different techniques vary with regard to sensitivity, selectivity and speed; furthermore measuring NH 3 is often expensive, extensive and time consuming [START_REF] Aneja | Summary of discussion and research re-commendations[END_REF][START_REF] Harper | Atmospheric ammonia: issues on transport and nitrogen isotope measurement[END_REF]. Nowadays among all the available techniques for measuring or estimating NH 3 , the most popular are (i) fluxes measurement approaches, as enclosure methods [START_REF] Mosier | Chamber and isotope techniques[END_REF] and micrometeorological methods [START_REF] Kaimal | Atmospheric Boundary Layer Flows: Their Structure and Measurement[END_REF], (ii) concentration-based dispersion modelling, Lagrangian [START_REF] Flesch | Deducing ground-toair emissions from observed trace gas concentrations: A field trial[END_REF], Eulerian [START_REF] Loubet | An inverse model to estimate ammonia emissions from fields[END_REF] or Gaussian [START_REF] Carozzi | A note on estimating the effect of a limited fetch on micrometeorological evaporation measurements[END_REF] types, and (iii) ammonia emission models (Gènermont and Cellier, 1997). The emission of NH 3 at field-scale depends on the interaction of various factors which contribute to decrease or increase the losses: fertiliser type (nitrogen content, pH, dry matter), soil type (water content, soil reaction), cultivation techniques (amount and application methods of fertilisers) and climatic conditions (temperature, wind speed, rainfall) [START_REF] Sommer | Effects of temperature, wind speed and air humidity on ammonia volatilization from surface applied cattle slurry[END_REF][START_REF] Moal | Ammonia volatilization following surface-applied pig and cattle slurry in France[END_REF][START_REF] Génermont | A mechanistic model for estimating ammonia volatilization from slurry applied to bare soil[END_REF]Sommer et al., 2001;Sommer and Hutchings 2001;[START_REF] Søgaard | Ammonia volatilization from fieldapplied animal slurry -the ALFAM model[END_REF][START_REF] Misselbrook | Dietary manipulation in dairy cattle: Laboratory experiments to assess the influence of ammonia emissions[END_REF]. The aim of this study was to estimate the NH 3 emissions caused by surface spreading of dairy slurry by using a concentration-based dispersion model. The quantification of the NH 3 fluxes was obtained by applying the backward Lagrangian Stochastic model (bLS) WindTrax [START_REF] Flesch | Backwardtime Lagrangian stochastic dispersion models, and their application to estimate gaseous emissions[END_REF]2004), since it has been increasingly employed in the last years (Sintermann et al. 2011b;[START_REF] Ni | The measurement time required for determining total NH 3 losses after field application of slurries by trail hoses[END_REF]. The model was implemented by the use of passive diffusion samplers [START_REF] Tang | Development and Types of Passive Samplers for Monitoring Atmospheric NO 2 and NH 3 Concentrations[END_REF] and the measure of atmospheric turbulence.

MATERIALS AND METHODS

The experimental site

The trial was performed from the 9 th to 17 th of October 2011 in Bigarello (MN), (Lat. 45°11' N, Long. 10°54' E, Alt. 23 m a.s.l.) in a field of 4.3 ha directions of the wind. The Monin-Obukhov length L [m] was derived from the Monin and Obukhov similarity theory (MOST, [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF] under horizontally homogeneous and steady state conditions:

[eq. 2]
where T [K] is the mean air temperature within the surface boundary layer, k is the von Kàrmàn's constant (0.41), g is the acceleration of gravity (9.81 m s -2 ) and w' T' is the covariance between w and T. The roughness length was derived from the wind speed profile relationship, as:

[eq. 3]
where U is the magnitude of the horizontal component of wind speed (m s -1 ) and Ψ m is a Monin-Obuhkhov universal function for momentum, estimated from the approach described by [START_REF] Flesch | Deducing ground-toair emissions from observed trace gas concentrations: A field trial[END_REF]. The u * and L were further filtered (u * > 0.2 m s -1 and |L| > 5 m) to guarantee the condition for the MOST application [START_REF] Flesch | Deducing ground-toair emissions from observed trace gas concentrations: A field trial[END_REF][START_REF] Hensen | Estimation of NH 3 emissions from a naturally ventilated livestock farm using local-scale atmospheric dispersion modelling[END_REF][START_REF] Loubet | Advection of NH 3 over a pasture field and its 32 Italian Journal of Agrometeorology -3[END_REF]. Moreover, to parameterise the bLS model, a constant value of z 0 (0.028 m) was calculated as the median value of eq. 3 over the experimental period.

The backward Lagrangian Stochastic model WindTrax

The WindTrax model (Thunder Beach Scientific, Halifax, Canada) is based on the backward Lagrangian stochastic dispersion theory described by [START_REF] Flesch | Backwardtime Lagrangian stochastic dispersion models, and their application to estimate gaseous emissions[END_REF]2004) and it has been employed to estimate the transfer coefficient D (s m -1 ). The transfer coefficient is used to derive the flux of NH 3 , S (μg m -2 s -1 ), emitted from the fertilised surface, on the basis of the NH 3 measured concentrations (C and C bgd , in μg m -3 ), from the relationship:

[eq. 4]
where D is retrieved by the model as the number of the interactions (N source ) between the source area and the thousands of trajectories (N)

The latter sampling point was used to measure the background level of NH 3 concentration (C bgd ). Samplers were exposed in three replicates at the height (z) of 1.25 m from the displacement height d = 0, corresponding to the ground. The positions of the ALPHA samplers and the shape of the fields were mapped using a GPS device. ALPHA were replaced a minimum of twice per day, after dawn and just before sunset, in order to monitor the change of atmospheric turbulence, which affects the dispersion of pollutants. During the daylight hours of the spreading day and the day after, the samplers replacement was done every three hours to have a more detailed time step. On the third day the replacement was done each 6 hours and subsequently every 12 hours. The exposed filters were leached with deionised water (3 mL) and then analysed by spectrometric detection (FIAstar 5000 system, FOSS, Denmark) through a gas semi permeable membrane (ISO 11732, 1997), in order to measure the concentration of NH 4 -N (mg L -1 ). The air NH 3 concentration (μg m -3 ) was then calculated by multiplying the NH 4 -N concentration, the volume of air sampled in one hour (V a = 0.003241315 m 3 h -1 ), the time of exposure (hours) (Sutton et al., 2001b;[START_REF] Tang | Analysis of temporal and spatial patterns of NH 3 and NH 4 + over the UK -2006[END_REF], and the stoichiometric ratio between NH 4 -N and NH 3 . The mean concentration and the standard deviation of the three replicates were calculated both for background (σ bgd ) and field measurements (σ C ).

Micrometeorological measurements

Micrometeorological measurements were performed to supply the parameters of atmospheric turbulence to the bLS model WindTrax (see section 2.4). The friction velocity (u * ), the Monin-Obukhov length (L) and the surface roughness length (z 0 ), together with wind speed (U) and wind direction (WD), were derived from a three-dimensional ultrasonic anemometer (USA-1, METEK GmbH, Elmshorn, Germany). The sampling frequency was 10 Hz and the device was set in the centre of the field at the same height of the ALPHA samplers. Friction velocity (m s -1 ) is derived from the Weber's formula (1999):

[eq. 1]
where u' and w' indicate the fluctuations of the wind components u and w along the three minimum values over the night time, where the amplitude and the magnitude of the peaks are directly related to the presence of horizontal wind. The negative peaks of z/L indicated the daily instability whereas positive peaks showed the nocturnal stability. Night and early mornings were characterized by phases of stability conditions, excepted for 14 th and 15 th of October where cloud cover and a strong and persistent wind occurred.

Ammonia concentrations

Fig. 2 shows the concentrations of NH 3 (μg m -3 ) measured over the experimental period in the centre of the field and present as background. Before the application of the fertiliser, the concentration was 1.7 μg m -3 higher than the background value. In the first three hours after the spreading, the field NH 3 concentration rose up to 505 μg m -3 , followed by a fast increase to 820 μg m - 3 in the subsequent three hours. In the following measurement period (6 to 15 h), the concentration decreased down to 300 μg m -3 and then increased during the night time hours to 539 μg m -3 . After the first 24 hours the concentration gradually decreased, reaching the background concentration after further 60 hours. In the last measurement period (84 to 168 h) the magnitude measured in the centre of the field and as background was coincident. Tab. 1 shows the statistics relative to such NH 3 concentrations. The variability (σ) of the three ammonia samplers ranged from 0.2 to 31.8 μg m -3 . In Fig. 3 the generated by the model and located upwind from the position of the two NH 3 samplers in the space (see eq. 5).

[eq. 5]

where (w 0 ) is the vertical wind speed of those trajectories that intersect the source area. The dispersion model used with the long term exposure samplers can be applied only with short time intervals (typically 30 min or 1 hour), because of the strong change of the atmospheric stability over a timescale of few hours. In the same way, the estimation of S is possible considering only the periods of atmospheric stationarity reached by short integration time of the turbulence parameters (u * , L and z 0 ), as the MOST theory states. For the determination procedure of S three hypotheses have to be assumed: (i) nonreactivity of the emitted NH 3 in the atmosphere, (ii) spatial homogeneity of the flux from the surface and (iii) steadiness of z 0 [START_REF] Loubet | Experimental assessment of atmospheric ammonia dispersion and short range dry deposition in a maize canopy[END_REF][START_REF] Loubet | Advection of NH 3 over a pasture field and its 32 Italian Journal of Agrometeorology -3[END_REF][START_REF] Nemitz | Aerosol fluxes and particle growth above managed grassland[END_REF].

RESULTS

Micrometeorological conditions

The trends of u * and the atmospheric stability parameter (z/L) measured together with their statistics, are shown in Fig. 1 and Tab. 1, respectively. Friction velocity marked the typical high peak values during the daylight hours and the

Ammonia emissions from slurry spreading

The flux estimated by the bLS model WindTrax and the error due to the dispersion of the NH 3 concentration mean value are displayed in Fig. 4. The emission trend showed high peaks immediately after the spreading, with a maximum value (163 μg m -2 s -1 ) reached after 3 hours and 20 minutes (see Tab. 1). Subsequently, the emission decreased quickly to rise again during the night time hours (15 to 21 h). The last high peak occurred after 24 hours because of the increase of the solar radiation. In the last measurement period (24 to 168 h) a gradual reduction to low values (2 μg m -2 s -1 ) was detected. Errors ranged from 0 to 4.7 μg m -2 s -1 , with a mean value of 0.4 μg m -2 s -1 .

Dynamics of soil NH 4 -N, NO 3 -N and pH

The dynamics of N into the soil as ammonium and nitrate form, together with the trends of SWC and pH at 0-5 and 5-15 cm depth, are displayed in Fig. 5a and5b, respectively. The measurement period ranged from -72 to 96 hours with a time step of 24 hours, where time 0 represented the time immediately before the manure spreading. Samplings at -48 and -24 hours were not carried out because of the high soil moisture. The content of NH 4 -N in the soil profile (0-15 cm) was constant before the slurry application (around 1 mg L -1 , from 0 to -72 hours) and rose up after the fertilisation, particularly in the first layer (0-5 cm). relation between σ and the mean values of concentration measured in the field and as background is displayed for each measurement period. The relation between C bkg and σ bgd was not significant (R 2 = 0.10, P>0.05), while a significant relationship was detected between C and its σ C (R 2 = 0.54, P<0.01). Subsequently, it quickly fell down because the concentration of TAN in soil surface decreased as consequence of emission itself, infiltration, absorption in the soil matrix, or nitrification (van der Molen et al., 1990;[START_REF] Sommer | Ammonia emission from mineral fertilizers and fertilized crops[END_REF]. Cumulative NH 3 loss observed at the end of the trial was visibly exhausted and equal to 30.2 kg N ha -1 , corresponding to an emission factor (EF) of 44.4% of the TAN applied. The emission ranged from 40% to 53% between 6 to 12 hours from the spreading. Similar findings were obtained by other authors, measuring a range of total loss between 30 and 70% of the TAN [START_REF] Sommer | Techniques and strategies for the reduction of ammonia emission from agriculture[END_REF][START_REF] Meisinger | Ammonia volatilization from dairy and poultry manure. Managing, Nutrients and Pathogens from Animal Agriculture[END_REF]. After 24 and 48 hours the total loss reached 73% and 88%, respectively, getting to the 100% of the emission at the end of the trial at 168 hours. The trend of the cumulate emission can be represented by the Michaelis-Menten equation (eq. 6) already used in this context by Søgaard et al. (2002): [eq. 6] where N max is the total loss of NH 3 (fraction of TAN applied) and K m [h] is the time t satisfying N(t) = ½ N max . A low value of K m indicates that a high proportion of the total NH 3 loss occurs soon after application, whereas a high K m value indicates that losses occur over a longer period. In our case after 10 hours the 50% of the total NH 3 was emitted (K m = 10 hours). The EF obtained in this field trial is coherent with what has been reported in the recent literature for similar conditions (see the review by Sintermann et al., 2011b). Furthermore, applying the regressive model ALFAM [START_REF] Søgaard | Ammonia volatilization from fieldapplied animal slurry -the ALFAM model[END_REF], based on the data deriving from 800 experiments and on the Michelis-Menten equation type described above, the EF obtained resulted 40% of the TAN. The fluxes estimated by the bLS model WindTrax showed high levels of emission when z/L parameter assumed negative values, or else when the atmosphere is in unstable condition. That was evident during the hours immediately close to the application of the slurry, characterized by high levels of concentration. The high fluxes observed during the night time hours of the spreading day (15 to 24 h) were due both by high levels of ammonia concentration (see Fig. 2) and the alternation of the atmosphere conditions. In fact, over the night, atmospheric condition can range from high stability, where the vertical gradients of concentration are enhanced to very small

The highest values in the two soil layers were observed 48 hours after the spreading. The values of NO 3 -N concentration did not seem to have a specific trend over time. In both layers the maximum values were reached at 72 h, showing the highest values in the lower layer (5-15 cm). At 72 and 96 hours the NO 3 -N concentration appeared to be higher in the 0-5 cm depth (+26% and +44%, respectively). The values of pH were reduced by 0.23 and 0.19 units in the 0-5 and 5-15 cm, respectively. The pH fast decrease started after 48 hours. The SWC showed values ranging from 24 to 31% in the 0-5 cm profile, with the maximum value observed 24 hours after the spreading. In the deeper layer a wider variation of SWC occurred (from 22% at -72 hours to 35% at 48 hours).

DISCUSSION

The NH 3 volatilisation estimated over the field trial followed the typical trend reported elsewhere (i.e. Sommer and Hutchings 2001;[START_REF] Powell | Dairy Slurry Application Method Impacts Ammonia Emission and Nitrate Leaching in No-Till Corn Silage[END_REF]. The highest emission rate was recorded immediately after the slurry application. 1983). Lastly the increase of the SWC due to the slurry application was observed in the first 24 hours in the 0-5 cm layer, whereas in the 5-15 cm layer the effect was extended till 48 hours.

CONCLUSIONS

Surface application of dairy slurry determined high ammonia emission levels within the first 24 hours, followed by a first rapid and then gradual decrease to low emission levels in the subsequent days. Ammonia losses were the 44% of the total TAN applied. The results obtained highlighted the need of low-emission techniques, such as surface spreading with incorporation and injection, which are supported by scientific data showing the reduction of ammonia emission under experimental conditions. However, the incorporation of slurry should be done as soon as possible because the 50% of the total ammonia emission occurred within 10 hours. Reliable NH 3 quantification has to be considered a central aspect for decision makers, to promote abatement strategies and to derive emission factors used in national and international emission inventories. diffusivity, to unstable conditions, in which concentration gradients are small due to the intense turbulent activity of the surface layer [START_REF] Erisman | Vertical distribution of gases and aerosols: the behavior of ammonia and related components in the lower atmosphere[END_REF][START_REF] Famulari | Development of a low-cost system for measuring conditional time-averaged gradients of SO 2 and NH 3[END_REF]. Another factor affecting ammonia emissions was the vertical wind speed (w), whereas the horizontal wind speed is not involved. In fact, as occurred from 14 th to 15 th October, the strong horizontal wind speed recorded did not affect the emissions that remained at low level. On the opposite, during the night time hours of the spreading day, when the value of vertical wind speed was high, an increase of emission occurred. The pattern of NH 3 emissions followed the trend of the measured concentrations and the high emission peaks took place in correspondence with the high values of concentration. Over the spreading daylight hours and the subsequent two days, the frequency of the ALPHA samplers change (3 hours in the first two days and 6 in the third, instead of 12 hours) allowed to obtain a more detailed ammonia concentrations and fluxes. In fact, integrating these concentrations over a time of 12 hours, the peaks of emission resulted lower and the final EF decreased by 5% (from 44 to 39% TAN). Therefore the temporal detail in the acquisition of the NH 3 concentration is a crucial point for a correct determination of the emission [START_REF] Carozzi | Dynamic of ammonia emission from urea spreading in Po Valley (Italy): relationship with nitrogen compounds in the soil[END_REF]. [START_REF] Loubet | An inverse model to estimate ammonia emissions from fields[END_REF] found that for an exposure time between 2 and 12 hours, the underestimation of NH 3 emission was not expected to be larger than 5 to 10% (± 5%) in a 100 m 2 surface area. Moreover, the use of inverse dispersion models, coupled with long term exposure concentration samplers, produced similar results if compared to the fastest and error-prone Eddy Covariance approach, when measuring NH 3 fluxes at field scale [START_REF] Ferrara | Eddy covariance measurement of ammonia fluxes: Comparison of high frequency correction methodologies[END_REF][START_REF] Carozzi | A note on estimating the effect of a limited fetch on micrometeorological evaporation measurements[END_REF]. Despite the low sampling frequency, an increase of the NH 4 -N concentration into the soil after the application of the slurry was recorded. As emissions ended also the NH 4 -N concentration decreased, whereas the NO 3 -N concentration started to increase due to the nitrification process (which was most evident at 0-5 cm depth soil profile). As reported by Gènermont and Cellier (1997) and [START_REF] Misselbrook | Dietary manipulation in dairy cattle: Laboratory experiments to assess the influence of ammonia emissions[END_REF] NH 3 emission and soil pH are directly related. In our case they started to decrease at the same time (after 24 h). Such behaviour was probably due to the release of protons (H + ) caused by the transition from NH 4 + to NH 3 in the nitrification process [START_REF] Van Breemen | Soil acidification from ammonium sulphate in forest canopy through fall[END_REF][START_REF] Freney | Volatilization of ammonia[END_REF]

Fig. 1

 1 Fig. 1 -a) friction velocity (u * ) and wind speed (grey dotted line) over the experiment; b) atmospheric stability parameter (z/L). Fig. 1 -a) velocità di frizione (u * ) e velocità del vento (linea grigia tratteggiata) per la durata dell'esperimento; b) parametro di stabilità atmosferica (z/L).

Fig. 4

 4 Fig. 4 -a) NH 3 flux simulated by the inverse dispersion model WindTrax. b) uncertainty in the modelling approaches employed due to the uncertainty in the concentrations measurements. Fig. 4 -a) flussi di NH 3 stimati attraverso l'impiego del modella a dispersione WindTrax. b) incertezza nella stima dei flussi attraverso il modello, dovuti alla incertezza nella misura delle concentrazioni.

Fig. 3 -

 3 Fig.3-Scatter plots reporting the average of NH 3 concentration and relative standard deviation (σ) for measurements in the field (∆) and at the background (▲). Fig.3-Grafico a dispersione riportante la media delle con centrazioni di NH 3 e la deviazione standard (σ) misurate nel centro del campo (∆) e come background (▲).

Fig. 3 -

 3 Fig.3-Scatter plots reporting the average of NH 3 concentration and relative standard deviation (σ) for measurements in the field (∆) and at the background (▲). Fig.3-Grafico a dispersione riportante la media delle con centrazioni di NH 3 e la deviazione standard (σ) misurate nel centro del campo (∆) e come background (▲).

Fig. 2 -

 2 Fig. 2 -NH 3 concentrations measured in the centre of the fertilised field (dark line) and in its background values (gray line). The error bars show standard deviations. Fig. 2 -Concentrazioni di NH 3 misurate nel centro del campo fertilizzato (linea nera) e come background (linea grigia). Le barre di errore mostrano le deviazioni standard.

Fig. 5 -

 5 Fig. 5 -Trend of soil N-NH 4 + , N-NO 3 -, pH and water content (SWC) during the time after slurry spreading in the field trial: a) 0 -5 cm depth; b) 5 -15 cm depth. Fig. 5 -Andamento del valore di pH e del contenuto di N-NH 4 + e N-NO 3 -e acqua (SWC) nel suolo durante le ore dello span dimento del liquame in campo: a) profilo 0 -5 cm di profondità; b) profilo 5-15 cm di profondità.

Italian Journal of Agrometeorology -3/2012Rivista Italiana di Agrometeorologia -3/2012