Root Loci Analysis of a Simple Quasi-Polynomial And Its Use For Relay-Based Parameter Identification of a Delayed Infinite-Dimensional Heat-Exchanger Process


 The focus of this contribution is twofold. The first part aims at the rigorous and complete analysis of pole loci of a simple delayed model, the characteristic function of which is represented by a quasi-polynomial with a non-delay and a delay parameter. The derived spectrum constitutes an infinite set, making it a suitable and simple-enough representative of even high-order process dynamics. The second part intends to apply the simple infinite-dimensional model for relay-based parameter identification of a more complex model of a heating-cooling process with heat exchangers. Processes of this type and construction are widely used in industry. The identification procedure has two substantial steps. The first one adopts the simple model with a low computational effort using the saturated relay that provides a more accurate estimation than the standard on/off test. Then, this result is transformed to the estimation of the initial characteristic equation parameters of the complex infinite-dimensional heat-exchanger model using the exact dominant-pole-loci assignment. The benefit of this technique is that multiple model parameters can be estimated under a single relay test. The second step attempts to estimate the remaining model parameters by various numerical optimization techniques and also to enhance all model parameters via the Autotune Variation Plus relay experiment for comparison. Although the obtained unordinary time and frequency domain responses may yield satisfactory results for control tasks, the identified model parameters may not reflect the actual values of process physical quantities.


Introduction
It is a well-known fact that dozens of industrial processes, including chemical ones, as well as social, economic, and other everyday systems are affected by latencies and delays [1][2][3][4][5][6] . Delays appear mainly due to mass, energy, and data transportation in the process and network interconnections, and their existence is closely related to distributed parameter systems. In modern discrete-time control systems, delays also arise from the human-machine interaction and signal sampling and processing 7 . As complex systems include internal feedback loops, internal delays must be considered along with the input-output ones; nevertheless, internal delays are often ignored when process modeling. However, such an approach can be unreasoning as the solution of partial differential equations (PDEs) -representing the reign of many industrial process models -often results in functions with lumped and distributed delays [8][9][10] .
On the other hand, time-delay models (TDMs) may be very good estimators of some systems and processes dynamics, even if any significant physical delay is not supposed to appear in the process. TDMs have the form of functional differential equations, or more specifically, delay differential equations (DDEs), instead of PDEs. Even a simple TDM can express the dynamics of a high-order non-delay model 11,12 with sufficient accuracy for control design. However, these models are infinite-dimensional because of the transcendental form of the characteristic equation (CE) 13 . All infinitely many solutions of the CE constitute the TDM spectrum of characteristic values (or poles). The pole loci most significantly determine the dynamic and stability properties of the model 14 . The infinite nature of the TDM spectrum yields its advantages and disadvantages when estimating the actual system dynamics. The dominant (i.e., usually the rightmost) subset of poles can match the pole loci of a high order system; however, one must be careful about other (uncontrolled) TDM poles, especially of highfrequency ones.
Various methods and techniques for the TDM spectrum analysis and its pole loci estimation exist; see a survey by Pekař and Gao 15 . We do let provide the reader with just a few. If TDMs have the so-called commensurate delays solely, pole loci can be determined analytically via the Lambert W function 16 . However, this is the only exact method besides the direct solution of the CE via the analysis of the distribution of the roots of the corresponding CE in the frequency domain, see, e.g., the work by Amrane et al. 17 . The family of numerical methods includes a wide variety of approaches and techniques that are based on, e.g., the mapping of real and imaginary parts of CEs solutions 18 , bifurcation analysis of DDEs 19 , full discretization of TDMs 3 , continuation property of TDMs pole loci 7 or on structural properties of a class of functional Vandermonde matrices 20 .
TDMs pole loci analysis is closely related to the inverse problem of the (partial) spectrum assignment, i.e., the determination of the model parameter (or even delay) values so that a subset of the characteristic values match the prescribed positions while other poles are sufficiently away from the chosen loci. There, again, exist several computational methods as well as principles of how to select the desired loci. The direct root assignment is essentially the most straightforward technique that gives rise to the solution of a set of algebraic (linear or nonlinear) equations 21 , the dimension of which is given by the number of assigned roots, their multiplicity and complexity. Its computational simplicity is, however, ransomed by a danger of a possible existence of model poles located right from the prescribed ones, which yields a problem of the poles' dominancy. Only a few tools ensure that the desired TDMs poles are dominant, e.g., a modified Nyquist stability criterion can be applied 22 , yet it is based on a graphical trial-andreset procedure. Recently, several ad-hoc results for single 17 and multiple real 20,23 prescribed poles or even single a complex conjugate pair 24 guarantying their dominancy have been derived; however, they usually levy large computational burden. Alternatively, the root dominance can be a posteriori checked using the argument principle (i.e., the Mikhailov curve based) approach 25 or via the solution of a special convolution integral 26 , which requires an advanced mathematical effort as well. Whenever the direct assignment is not satisfactory, a numerical spectrum optimization can be made, e.g., by the quasi-continuous shifting of the roots 27,28 or using its combination with the minimization of a specific fitness function reflecting the remaining spectrum, robustness issues, etc. [29][30][31] . Unfortunately, a non-convex optimization problem must be solved in many cases, see, e.g., 32 and references therein.
The use of relay in the feedback control system represents a favorite system parameters' identification and automatic controller tuning framework that have received a great deal of attention since the pioneering work by Åstrom and Hägglund 33 , where the ideal on/off relay was used to generate sustained (ultimate) oscillations. This parameter estimation framework enables to prevent the process output from drifting too far away from the reference signal (setpoint), which is required for many industrial processes. During recent decades, a multitude of derived techniques and methods have been developed 34,35 that have found a great favor of practitioners, especially in chemical and process engineering [36][37][38][39][40][41] .
Three families of approaches to evaluate unknown model parameters 34 exist. Namely, using a describing function (DF) represents the most common approach 33,38,42 . Roughly speaking, this function is a linear approximation (usually based on the Fourier series expansion) of the nonlinear relay behavior. Second, the curve fitting approach attempts to fit the feedback response in the time domain based on an analytic formulation of the response [43][44][45][46] . As third, the frequency fitting does the same yet in the frequency domain, which corresponds to the seeking of multiple frequency points, besides the ultimate case [47][48][49] .
It is known that a sufficiently accurate process model can reduce errors in controller tuning 50 . The ultimate gains obtained from the standard (on/off) relay feedback oscillation amplitudes can have errors of over 15% 51 . Besides, Jeon et al. 52 pointed out that model parameters obtained from the sustained relay oscillations can be insufficient if there is a mismatch in the model order and process dynamics, which gives rise to the need for more complex models. However, as the original method suffers from significant errors in model parameters' estimation and only one (critical) point of the frequency characteristics (i.e., two model parameters) can be obtained from the test, researchers have developed methods to fit the parameters more precisely and/or to search for more frequency response points under one or more relay tests.
Regarding the former group of methods, improved accuracy can be obtained by compensating for the phase lag caused by the relay module 53 -which is suitable for higherorder processes or those with large input-output delay, by using a biased relay 54 , a relay with two-band hysteresis to reduce the oscillation frequencies 55 or relays with multiple switching 56,57 that prune the relay oscillation harmonics and the effect of the input nonlinearity. Unfortunately, the reduction of the oscillation frequencies increases the experimental times. Other techniques, e.g., attempt to obtain as sinus-like relay output as possible by using a relay with saturation 35,58,59 , to reduce the effect of noises and disturbances 50,60 , to apply the so-called area methods that integrate specifically modified time responses 61 or to use asymmetrical limit cycle 62 or a shape factor 63 . However, most methods suffer from the sensitivity to plant-model structural mismatch 58 .
The simplest approach how to obtain more frequency points is to perform more than one relay test using an additional integrator 64,65 , via the parasitic relay 66 , a relay with hysteresis 67 or the biased relay 42 that alter the DF 68 . Li et al. 69 proposed a well-applicable method called the Autotune Variation (ATV) that introduces an artificial delay in an additional test following the standard relay experiment. The technique was further improved by Kim 70 , Marchetti and Scali 38 , and Scali et al. 71 . These methods are popular among practitioners due to their computational simplicity. However, multiple experiments may be time consumptive; therefore, researchers have attempted to gain information about more frequency points or other system dynamic features under a single relay test. Even a purposefully induced asymmetry can be used to determine additional frequency response points without performing additional relay feedback tests 72 . Nevertheless, the given asymmetry may yield the termination of oscillations, which fails the relay test. Numerous other techniques exist, such as the shifting method 51,73,74 , the use of relay transient 48,54,75 or the computation of weighting moments 47,61 .
Most of the linear process models used in literature when applying some of the relay-based identification experiments include input-output delays within First-Order plus Dead Time (FOPDT) or Second-Order plus Dead Time SOPDT models, either stable 44,45,46,50,55,63,69,76 or unstable 43,44,51,[77][78][79][80][81] ones. The models are represented simply by a series of a delay-free loworder submodel and the delay element. However, FOPDT models can suffer from poor performances for some low-order processes with fast parasitic dynamics 61 . Besides, the delay value is hardly ever considered as significant or dominant, even if such an assumption is often unreasonable in practice 82 .
Surprisingly, not many results concern internal delays incorporated in TDMs. Parameter estimate of TDMs is more complicated than parameter estimate of standard models 47 . The pioneering work in this field was presented by Vyhlídal and Zítek 12 . The authors comprised a first-order-derivative model with one internal delay parameter within the standard relay test 33 .
We do let call this TDM the Simple First-Order one (SFOTDM) for further presentation. The identified SFOTDM was further used for internal model control design. These results were taken as a starting point for parameter estimation of a TDM of the continuous stirred-tank reactor. The parameter values are further enhanced by solving nonlinear objective functions governed by the difference of the model and measured responses in the time domain via anoptimization algorithm. Nevertheless, a relay was not used in the proposed design; the authors only referred to the awkward need to determine the input-output delay value before the relay test. Pekař and Prokop 83 compared the use of the saturated relay 42 and the limit-cycle evaluation using the exponential decaying followed by the discrete-time Fourier transform 54,66 . The authors considered a first-order-derivative TDM that included three non-delay parameters, one internal delay parameter, and the dead time as the model of a circuit system with heat exchangers 84 . An artificial delay was used for the additional relay test. In 47 , the method of moments 64 was applied to the serial combination of the SFOTDM and a low-pass filter represented by the delay-free first-order submodel 85 , and to its high-order generalization that, however, does not represent a universal TDM. Two versions of the shifting technique (or, the shift transformation) were applied to the preceding model by Hofreiter 74,86 . The DF was based on the fundamental-harmonic (i.e., higher harmonics are neglected) Fourier series expansion of the shifted control (input) signal. A biased relay with hysteresis was used because of practical reasons; however, it was not explicitly included in the algorithm.
As this study concerns a circuit process with heat exchangers, it should also be noted that relay-based identification experiments (followed by a controller design in many cases) have been applied to heat-exchanger systems. For instance, a case study on an autotuning control method for a cross-flow heat exchanger was published in 87 . Jin et al. 88 presented a Ziegler-Nichols-based method based on using the ultimate gain (instead of on a nonlinear element) to get the sustained oscillations. Some researchers use nonlinear models such as Wiener-type or Hammerstein-type 89 . The reader is also referred to 37 and references therein.
Let us provide the reader with the motivation to perform the presented research. Models of industrial processes usually include a large number of unknown parameters. Hence, when applying the above-celebrated relay-based identification tests in practice, multiple experiments and/or solutions of nonlinear optimization problems are usually needed (even for linear models). As mentioned above, we focus on a circuit process with heat exchangers with large input-output and internal delays 84 . Considering the simplest single-input single-output case, the derived model includes six non-delay parameters, one internal delay parameter, and two parameters in the input-output relation. Let us call the model as Heat-Exchanger TDM (HETDM) for further presentation. Therefore, any attempt to fit these nine parameter values (e.g., in the frequency domain) requires a good initial guess. Hence, the main idea of the presented research is to perform a two-step relay-feedback identification procedure. The SFOTDM is assumed in the first step, which requires a relatively low computational effort when estimating model parameters' values. The second step adopts the identified SFOTDM in the sense that dominant characteristic values (poles) of the model coincide with the dominant poles of the analytically-derived HETDM. The pole assignment yields the determination of the characteristic function parameters. These values can either be fixed while remaining model parameters are then set from the single experiment data or they constitute the initial estimation that are further enhanced along with other undetermined parameters that do not depend on poles. As the bridge between the two basic steps, a thorough analysis of pole loci of the SFOTDM resulting in explicit and implicit formulae and a simple graphical procedure is made.
This quasi-polynomial root analysis constitutes a substantial contribution to the presented paper.
Relay feedback experiments for both the used models use the standard (i.e., on/off) biased relay for the initial estimation of the ultimate gain and the computation of the process static gain 90 . However, the selected asymmetry is small enough not to affect the applied DF yet sufficient to determine the static gain. The relay with saturation 42 is applied to enhance the DF evaluation. This nonlinear element can estimate the DF precisely in the ideal case. In the second step (i.e., for the HETDM), single and multiple relay tests are made. The single test is performed to identify the transfer function numerator parameters that are not affected by pole loci.
Contrariwise, the multiple test attempts to determine four frequency points via three additional experiments utilizing an artificial delay as per the ATV+ technique 38,69,71 . The response ultimate data and the given DF are processed via the well-established Levenberg-Marquardt (LM) method 91 and the Nelder-Mead (NM) algorithm 92 to solve a certain nonlinear frequency-based constrained optimization problem. Several numerical scenarios are benchmarked. Namely, LM and NM techniques are compared when using a single relay test to determine transfer function numerator parameters. In addition, the NM algorithm is used to estimate eight model parameters when applying the ATV+ test.
Contributions of the presented research can be summarized as follows: 1) Exact analytic rules to determine pole loci of the SFOTDM are derived.
2) Saturated relay feedback experiment is performed on the SFOTDM. The detected pole loci are then set to the initial parameters' guess for the HETDM (as a sufficiently accurate mathematical model of the circuit system with heat exchangers), which is enhanced via the LM method, under the single relay-test data.
3) Three scenarios are compared to determine the remaining model parameters and further enhance the already estimated ones via the ATV+ technique and the solution of a nonlinear optimization problem using the LM and NM algorithms. 4) An independent determination of numerator and denominator transfer function coefficients along with the pole loci assignment enables to reduce the number of necessary relay test for some of the scenarios.
The rest of the paper is organized as follows. "Methods and techniques" summarizes theoretical fundamentals of (retarded) TDM spectrum and feedback relay-based experiment, the model parameters' identification using a DF, and the LM and NM methods. "Results" has two fundamental subsections. The first subsection provides the reader with a detailed analysis of the SFOTDM pole loci. The second one presents the HETDM and all steps of determining its parameter values. Namely, the mathematical model of the HETDM is introduced, then the reader is acquainted with the poles assignment, the transfer function numerator estimation using a single relay test, and the complete model parameters estimation via the ATV+ technique. In "Discussion", the obtained results are discussed. Finally, "Conclusions" concludes the paper.
The standard notation is used throughout the paper, i.e., ,, denote the sets of complex, natural (excluding zero) and real numbers, respectively, n  expresses the n-dimensional

Relay-based parameter identification
As introduced above, experimental plant identification using the relay (or another simple      , and 0   expresses the hysteresis parameter. The , see, e.g. 67,68,73 . Then, the ideal on/off relay gives rise to In practice, the setting 0   is suitable when the feedback signal is affected by noise so that the switching relay rate can be reduced. The advantage of the option 0   lies, i.a., in the possibility to estimate the process static gain     for 0 t satisfying that sustained oscillations start for some 0 tt  35 .
Purposefully induced asymmetry can also be used to estimate and attenuate the load disturbance 42 . However, it may stop oscillations so that the relay test fails. In addition, model parameter identification with asymmetric relay yields an estimation error of up to 40% in a FOPDT case 50 .

Relay with saturation
Estimating the critical point at u  (or any other osc  ) does not provide an accurate enough parameter estimation for some processes, e.g., for those with significant time delays. For instance, an error of 23% for FOPDT models was reported 59 .
Model parameters identification can be improved by using saturation relay 35 Ideally, if the gain sat k is set optimally (i.e., AA  ), input and output signals has the same shape; hence, the DF

ATV+ technique
The disadvantage of this technique is the prolongation of the relay feedback experiment.
However, if the initial conditions are sustained oscillations, it lasts a significantly shorter time to restore the oscillations than starting from a constant steady state.

Parameter optimization methods
To solve (2) for given roots k s , (4), and (8), two well-established optimization algorithms are adopted. Their concise description to acquaint the reader follows.

Levenberg-Marquardt method
Consider a set of nonlinear differentiable functions where   0, 0, , 0 J means the Jacobian of f with respect to p , k expresses the iteration step, and 0   is an adjustable parameter (the so-called damping factor) 93 . Solution (11) The value of  (the so-called damping factor) may vary during iterations. One of the framework strategies is to decrease its value as the residual sum on the right-hand side of (13) (

 
Res p ) decreases, and vice versa. Let us introduce the multiplicative factor 0 . Particular choices of 11 ,  p , and  are discussed in "Relay-based parameter identification of heat-exchanger process".
A disadvantage of the LM algorithm is that the solution may converge to a local minimum (as for other Newton-type methods) or it may even diverge (especially, if  is set inappropriately).

Nelder-Mead method
Assume an unconstrained optimization problem, first The idea of the NM method 92 is to iteratively search for the optimal solution by moving a variable-shape simplex in the space of p . The simplex vertices represent the so-called test points. Once the initial simplex   The worst-valued vertex is reflected through k c p as Then, four scenarios can happen: shrink the simplex using (16).
Then 1 kk , re-order simplex vertices, and calculate (15), etc. If, however, inequality constraints on a subset  pp % are required, one may use the concept of barrier functions. That is, instead of the objective function   f p as in (14), the extended function    p is subject to the optimization procedure where 0 must be sufficiently small as soon as all

Root loci analysis of the simple quasi-polynomial
In this subsection, a thorough zero loci analysis of the SFOTDM 12 , and they determined ranges in which the model is asymptotically stable, aperiodic, and periodic. Moreover, intersections of pole loci trajectories with the imaginary axis for the generalized model were determined. Analogous conditions for which the model is stable, overdamped, critically damped, and underdamped were presented in 95 .
Asymptotic behavior of pole loci trajectories in infinity and nearby the imaginary axis were also studied in 96 .
Hence, our aim is to analyze the solution (or its rightmost subset) of the CE in and provide the reader with a simple guide how to compute these pole loci.
The first formula of (26) can be rewritten as The latter condition in (26) agrees with By comparison (27) and (28), it can be deduced that 1   is the only double real root. Now we must show that it is the only multiple root for any finite  . Such a root must satisfy (25) and simultaneously Since 0   , only two solutions of (29) exist The former solution in (30) yields the stability border due to Lemma 1 and   from (27) and (28), which indicates a root at infinity. It is, however, in contradiction with (29). The latter possibility in (30) gives    from (27) but    from (28), which yields a contradiction again.
As an alternative of the proof, one can easily deduce that (28) and (29) Fig. 6. ■

□
Proof. Regarding item a), the given range of  yields the rightmost complex solutions of (20) according to Theorem 1. Besides, its abscissa is within the range (24). All complex roots have their loci given by (31) and (32)    Despite its simplicity, the mathematical formulation of the HETDM and especially its dynamic properties are remarkable due to the model transcendental characteristic equation 84 .
As the model is multivariable, the relation between the heater power input   ut and the cooler outlet heat fluid temperature   yt is selected as the most interesting input-output pair. Note that both quantities are considered as their deviations from a steady state. The analytically modelled linearized relation reads

y t a y t a y t a y t b u t b u t
which is a DDE where 0 0 2 1 0 0 , , , , , b b a a a a Let us take these data as a benchmark for the significantly more straightforward relay-based experiment. As the values arise from determined physical quantities of the process, they are assumed to be closed to the actual (true) real-life values.

Simple model parameter estimation using the relay-based experiment
The first step of the identification chain is estimating the SFOTDM parameters, especially those of SFOTDM q (20). It has three substeps. First, the on/off relay with 0   (and  being sufficiently small), see (5), is used to estimate the static gain k in (19) as per (6). Second, the ideal relay ( 0   ) is applied to get the initial estimation of oscillation data and the input/output delay value. Finally, the saturation relay is used to improve the accuracy of oscillation parameters, which yields the SFOTDM parameters from (4) and (7). All the substeps can be done within a single experiment, saving time since the transition from particular substantial oscillations to others takes less time than setting the oscillations from a constant steady state. Let us use the notation , ss      for (19) to distinguish the SFOTDM parameters from those of the HETDM (for which no subscript is used). The combination of (4) and (7) (4) and (7) have to be used instead of (47).
, and 5 10    . The relay-test responses are displayed in Fig. 9. The arrows indicate when a particular relay starts to be used.   Fig. 10. The models in Table 2 are also equipped with Roman numerals to label them.   Note that the terminal frequency value in Fig. 11 Table 3 prove the above-introduced assumption that Models V and VI estimate the process frequency response reasonably on low frequencies (i.e., up to the oscillation frequency); however, they fail for higher ones. It is also worth noting that the estimation of the ultimate frequency is quite accurate since the true one is about 2 1.67 10

Parameter estimation of the heat-exchanger process model via pole assignment
Now, models (results) in Table 2 serve as the initial estimates for HETDM parameters identification. Several scenarios are tested and compared after the HETDM dominant spectrum assignment according to the pole loci of the SFOTDM.

Initial denominator parameters estimation
The rightmost spectra of SFOTDMs (see Table 2) poles are displayed in , , , , s s s s s . However, it is not always possible. For instance, SFOTDM I has only two significant real poles (the rest are too far in the LHP), and complex conjugate pairs of other spectra cannot be decoupled. Hence, only less number than five of SFOTDM poles are taken in (50). It is also worth noting that whenever \ k s  , the HETDM q in (50) is split into its real and imaginary parts.
We use the LM method (see "Levenberg-Marquardt method") for solving the pole assignment problem. Let us discuss the ideal initial parameters' selection,    (13)) are provided to the reader in Table 5. Corresponding dominant pole loci of the models are enumerated in Table 6.
Note that the QPmR software package 18 is utilized to compute the poles here.     Table 6. Dominant pole HETDM spectra.
Now three different scenarios of how to set numerator parameters of the model transfer function (45) or even eventually alter all the transfer function parameters based on relayexperiment data follow.

Numerator parameters estimation using Levenberg-Marquardt method
The first scenario adopts the LM method and utilizes data solely from the single relay test (i.e., without a necessity to perform additional experiments) that should indicate the ultimate data (i.e., the critical point of the Nyquist curve). Once the parameters are determined as in "Initial denominator parameters estimation", they are fixed. Hence, the HETDM transfer function numerator parameters are estimated to comply with conditions (4). The advantage of this scenario lies in fact that all the model parameters are found within a single test.
In more detail, the set of nonlinear algebraic equations to be solved reads where osc  holds when using the relay with saturation. As only two parameters can be determined by solving (53), two other ones need to be set a priori. Let s   be fixed as in Table 2. As the static gain is known, the following equality is substituted to (53)   Hence, the parameter set in (53) eventually reads   00 , and 0 b is then calculated using (54). The LM control parameters are set as in "Initial denominator parameters estimation".  Table 7).
Compared to SFOTDMs (Tables 2 and 3 Table 8 are very close to those displayed in Table 3.
The unit step responses and Nyquist plots of some selected HETDMs are displayed in Figs where the HETDM characteristic RQPs are fixed as in Table 5  Gs is used in (55) rather than that with the amplitude and phase since numerical experiments give better results (in the sense of (4)).
The most outstanding results are provided in Table 9 (two parameters' sets from each of the six pole spectrum families are selected again). Table 10 displays corresponding performance measures and implies that the accuracy of HETDMs in Table 9 is very close (or slightly worse) to models obtained in "Numerator parameters estimation using Levenberg-Marquardt method" (except for models from family I). This means that the models give better accuracy than SFOTDMs in the time domain yet only comparable ones in the frequency domain.
Unit step responses and Nyquist plots of selected models are given in Figures 14 and 15, respectively. Note again that other characteristics are close to some displayed ones. In the time domain, responses for model families I, II, and III almost coincide when I-a-2 is the fastest and II-c-2 is the slowest. The same assertion holds for families V and VI (VI-b-2 and V-a-2 give the fastest and the slowest response, respectively), while models IV significantly differ from the others. Nyquist plots of III-c-2 and III-f-2 are closest to each other at low frequencies and the pair II-c-2/III-f-2 at high ones. The pair VI-b-2/VI-d-2 almost coincides at the whole frequency range. Notice that although perfect optimization based on the measured ultimate data (see Table 10) is reached, model IV-a-2 does not provide excellent frequency response.   However, the curve for the latter model is closer to the original Nyquist plot than that for the former one. This implies that the model accuracy cannot be judged solely on the shape of characteristics.
Another question is whether optimizing more than three transfer function numerator parameters can improve the HETDMs accuracy.

Numerator/denominator parameters estimation using Nelder-Mead method via Autotune Variation Plus experiment
By substituting (54) into (45) a a a a  f b b  a a a a   f b b  a a Table 11).
Apparently, HETDM parameter identification based on the estimation of four Nyquist plot points results in significantly improved time-domain responses compared to 3-parameter optimizations (see "Numerator parameters estimation using Levenberg-Marquardt method" and "Numerator parameters estimation using Nelder-Mead method from single relay test").
Regarding the model accuracy in the frequency domain, substantial enhancement is achieved for low frequencies, which, however, does not hold for the whole frequency range. As can be deduced from Tables 8, 10,  In Figure 17, nearly nonsmooth step response of HETDM I-a-3 can be observed. An oscillatory response with high-frequency modes of nonnegligible amplitude can be seen for HETDM 4-a-3A.
The remaining (non-displayed) dynamic responses have a similar shape to the displayed ones; however, the plots are not as close as the 3-parameter optimizations. An exception appears for III-f-3, V-c-3, and VI-a-3 where nonminimum-phase-like time responses appear, yet the corresponding Nyquist plots do not prove this feature. The significant step response differences come from diverse input-output delays. Regarding Nyquist plots, the closest curves can be observed for pairs I-a-3/I-b-3 and IV-a-3A/V-b-3 at low frequencies. At higher frequencies, responses differ more significantly, especially those for HETDMs III-c-3, III-f-3, and V-c-3 are far from the remaining ones.

Discussion
Let us discuss observations made during the entire HETDM identification procedure and also point out some practical issues.
In our experiments, we have supposed that the hysteresis of the on/off relay is negligible, 0   . However, a suitable nonzero value has to be set in practice due to the measurement noise.
Such a setting prevents the relay from switching too frequently, which may cause failures.
Another important practical issue is the static gain estimation, according to (6). Whenever an asymmetry ( BB   ) is induced, the output of the feedback system also becomes asymmetrical.
The problem is that the original output setpoint then shifts, which may cause an erroneous estimation of the static gain. The setpoint shifting can be caused by the feedback nature of a relay experiment, process nonlinearities, and/or disturbances. If someone is unsure about the static gain, the step response test can be made. It is worth noting that disturbances also induce asymmetry of the ideal on/off relay experiment. In such a case, various methods of restoring symmetry can be applied 50 .
As a DF generally represents a linear approximation of a nonlinear element, it is impossible to estimate the critical (or another frequency) point exactly by nature. More precisely, neither the found frequency osc  nor the loci  meets the particular values of the actual (measured) process frequency response. This implies that even if the solution of (4) is perfect (see, e.g., HETDM V-a-2 in Table 10), the model does not provide sufficient results from the identification point of view. Besides, even if one or more points of the Nyquist curve are estimated well, the remaining course of the plot may vary from the desired loci significantly (see, e.g., HETDM I-a-3 in Fig. 18).
As can be seen from Tables 7, 9, and 11, relatively high ratios 00 / bb  , 00 / aa  , 00 / aa  ,    often occur. This unpleasant feature yields erroneous steady state computation or numerical instability when simulations (i.e., solving differential equations) due to digital representation of values in computer. Therefore, only some eventual models can be used for control system design and its verification.
Displayed step responses indicate that the initial input/output delay estimation 136.7 s   is pretty good. It can be observed that corresponding model families I, II, and III (see Tables 2, 3 The eventual SFOTDMs identified using relay feedback experiment were proved to be sufficient for controller design 12 . By matching SFOTDMs dominant pole loci with those of HETDMs and calculating remaining model parameters, performance measures very close to those of SFOTDMs have been obtained in this study. This implies that the eventual HETDMs based on the data from the single relay test (i.e., without artificial delays) can also be used for control tasks. However, the models seem to be insufficient from the identification point of view.
Fortunately, the use of the ATV+ experiment has brought about much improved models. As many diverse results have been computed with more or less the same cost function values, the identification problem for HETDMS seems to be multimodal task. Therefore, none of the eventual models (see Table 11) approaches the true parameter values (46) given by physical analysis of the thermal process.
The relay-based experiment can be improved in several ways. For instance, one has to be more careful when setting sat k of the saturated relay (see Fig. 4). The ideal setting should result in sinusoidal-like sustained oscillations of   ut. However, Figs. 9 and 16 indicate that sat k was set too high. On the other hand, one has to be aware of the necessary existence of sustained oscillations. Another way how to enhance the coefficient value estimation is to capture multiple points of the Nyquist plot, which may yield better matching of process and model curves for a wide frequency range [47][48][49] .

Conclusions
This study should have examined whether it is reasonable to identify parameters of a complex model of a thermal circuit system with internal delays by a parameter identification of a simpler delayed model followed by the models' poles matching. As the identification tool, the standard on/off relay with biased and unbiased feedback test and the relay with saturation have been used. The latter relay should have yielded a more accurate estimation of points on the frequency curve corresponding to sustained oscillations data. Once the simple model is found under a single feedback experiment, its dominant pole loci (of an infinite spectrum) are matched to those of the complex model, giving rise to the characteristic quasi-polynomial coefficients. A simple graphical-based method has been analytically derived to find these loci.
The Levenberg-Marquardt method has been applied to solve the pole assignment task.
Surprisingly, although only a few poles have been prescribed, some other uncontrolled poles have also been matched. Based on the single-test data, the remaining model parameters have been estimated by the solution of a nonlinear optimization problem (using the Nelder-Mead vs. the Levenberg-Marquardt methods). It has been proved that both the models have had similar time and frequency domain performances. While the eventual models may be sufficient from the control point of view, they fail regarding the accuracy of identified parameters. On the other hand, the proposed procedure enables the estimation of multiple parameters under a single relay test, which is its main benefit. However, we have also performed the ATV+ test with artificial delays to get multiple relay-feedback data, which has resulted in much better eventual models.
In the Discussion section, we have touched on some issues that have to be considered in practice and proposed possible further improvements to the proposed concept. Besides, one may apply advanced and more sophisticated optimization approaches to solve the tasks raised in this study. For instance, metaheuristic methods can be benchmarked in the future 98 . In addition, real-life experiments will be made to prove the concept.

Data availability
Data are available by L.P. upon request.
The authors declare no competing interests. Figure 1. A framework scheme of the relay feedback identification experiment.                     Table 4. Dominant pole SFOTDMs spectra. Table 5. HETDM characteristic RQP parameters after pole assignment. Table 6. Dominant pole HETDM spectra. Table 7. HETDM transfer function numerator parameters computed via the LM method.  Table 7). Table 9. HETDM transfer function numerator parameters computed via the NM method. Table 10. HETDM time-domain and frequency-domain model errors (ad Table 9).   Table 11).