Sampling Rate Optimization and Execution Time Analysis for Two-Degrees-of-Freedom Control Systems ^†

Abstract

The current journal paper proposes an end-to-end analysis for the numerical implementation of a two-degrees-of-freedom (2DOF) control structure, starting from the sampling rate selection mechanism via a quasi-optimal manner, along with the estimation of the worst-case execution time (WCET) for the specified controller. For the sampling rate selection, the classical Shannon–Nyquist sampling theorem is replaced by an optimization problem that encompasses the trade-off between the fidelity of the controllers’ representation, along with the fidelity of the resulting closed-loop systems, and the implementation difficulty of the controllers. Additionally, the WCET analysis can be seen as a verification step before automatic code generation, a computational model being provided. The proposed computational model encompasses infinite-impulse response (IIR) and finite-impulse response (FIR) filter models for the controller implementation, along with additional relevant phenomena being discussed, such as saturation, signal scaling and anti-windup techniques. All proposed results will be illustrated on a DC motor benchmark control problem.

1. Introduction

1.1. Literature Review

The problem of sampling rate choice to numerically implement a controller designed in the continuous-time domain has significant importance. A sub-optimal value of the sampling period represents a trade-off between the fidelity of the response, given by a shorter sampling rate, and implementability, obtained using a larger sampling rate. The most common manner to choose this sampling period is given by the Shannon sampling theorem [1]. However, this can only be a starting point in a control context, because, in practice, it can lead to unacceptable behavior of the discrete-time controller compared to its continuous-time equivalent, a significant justification being the presence of quantization errors unaccounted for in the signal and system models. In the available literature, the problem of choosing the sampling rate is specifically formulated for the purpose of the particular control system at hand. For example, for a structure with a P-type regulator proposed in [2], the key points used in sampling rate selection are the overshoot and the rise time. For the case of a PID-based control structure, an optimization criterion has been proposed in [3] to find the optimal and/or sub-optimal value of the sampling rate. The problem of sampling rate selection for parameter estimation of an induction machine has been solved via a metaheuristic procedure in [4]. A similar thematic is addressed in [5] in the case of permanent magnet synchronous motors using field programmable gate array devices, while the rapid control prototyping principle is illustrated for discrete-time closed-loop control of brushless DC motors in [6,7], respectively. Moreover, an important result regarding the dependence between the stability and the performance of the closed-loop system, and the choice of sampling rate is presented in [8].

After the first step of selecting the sampling period, an analysis regarding the practical implementation of the proposed control structure on a microprocessor-based system should be performed. This analysis needs to encompass several aspects which are not explicitly modeled in the control law, such as the register word lengths of the coefficients which directly influence the number of necessary assembly instructions, and saturation, overflow, underflow verifications on the computed command signals. In the specific context of rapid control prototyping (RCP), one key step before the proper code generation is to certify if its execution abides by the time span given by the sampling period. Moreover, other important verifications for phenomena such as underflow and overflow, or saturation and anti-windup techniques should be also performed before the code-generation step. As such, one important goal consists of approximating the number of software operations necessary to implement the designed control law, mandatory in the context of hard real-time systems. A survey of worst-case execution time computational methods is presented in the seminal paper [9]. A simplified study applied for control system-related interrupt service routines for state-space controller representations is presented in [10]. Another possibility for the computational error analysis is presented in [11].

In the Control Systems domain, the linear control system field plays an important role. Linear controller design techniques are available for both linear plants and nonlinear process models alike. Starting from the classical proportional–integral–derivative (PID) controller [12] towards the domain of robust controller [13,14,15], these techniques can be extended for nonlinear systems via the gain-scheduling method [16] or an additional passivity-based component [17]. Moreover, with great advantage in also coping with disturbance rejection problems, alongside reference tracking specifications are the so-called two-degrees-of-freedom (2DOF) control schemes [13]. In ref. [18], a 2DOF PID controller has been proposed, where the serial compensator is a classical PID controller, while the feedforward compensator is a PD controller, the integral effect being excluded due to the stability requirements. For the case of unstable plants with time delays, a discrete-time 2DOF control scheme has been proposed in [19]. As such, for reference tracking, the serial controller was designed using the ${\mathcal{H}}_2$ optimal control framework, while the feedforward controller has been tuned by imposing the desired closed-loop transfer function. Additionally, in terms of software toolbox implementations, the robust advanced PID (RaPID) toolbox described in [20] presents a set of possibilities to design a 2DOF PID control structure by minimizing certain error criteria, such as the integral of absolute error or integral of time multiplied by the absolute value of error.

1.2. Contributions and Paper Structure

The current journal paper represents an extension of the conference papers [10,21], and proposes a joint analysis on sampling time selection and execution time framing of the microcontroller operations into the previously established optimum for the case of a 2DOF control structure. The main contributions of the paper are:

(i)

to extend the functionals initially proposed in [21] in order to encompass the fidelity of both controllers from the 2DOF structure and of the resulting closed-loop systems, i.e., ensuring the tracking and servo behavior of the initially designed continuous-time controllers, along with the implementation difficulty of the proposed controllers in terms of quantization and sampling rate span;

(ii)

to formulate the optimization problem for sampling rate selection of a 2DOF control structure without considering any further normalizations of the final functional’s terms and propose a metaheuristic-based solution to find an optimal or quasi-optimal value of the sampling rate which can offer a good trade-off between implementability and fidelity;

(iii)

to extend the WCET analysis, initially proposed in [10] for the case of state-space regulator implementations, for the case of controllers implemented as infinite-impulse response (IIR) or finite-impulse response (FIR) filters, offering an upper bound estimation of the time span necessary to perform all the operations involved in the numerical implementation of the proposed 2DOF controller;

(iv)

to offer an exhaustive analysis which can be further performed in an automatic manner and to be integrated in our open-source MATLAB-based CACSD toolbox [22];

(v)

to present an end-to-end design procedure for the case of a DC motor control problem, starting from the controller design, followed by choice of the quasi-optimal sampling rate and the worst-case execution time analysis, finalizing with a detailed discussion.

The rest of the paper is organized as follows: Section 2 presents the 2DOF structure which will be studied, along with a mathematical background for the sampling rate choice and worst-case execution time analysis, and a set of theoretical results and optimization problems extended and adapted for the proposed control structure; in Section 3 a case study consisting of an end-to-end design procedure of a 2DOF PID controller is illustrated, while Section 4 deals with a set of conclusions and further research directions.

1.3. Notations

We denote by $\mathcal{C}\left(z,r\right)$ the circle with the center in $z\in \mathbb{C}$ and radius $r>0$. Additionally, $\mathcal{P}\left(G\right)$ and $\mathcal{Z}\left(G\right)$ denote the pole and zero sets of an LTI system G, respectively. An arbitrary sampling period will be denoted $T_s\equiv T>0$ throughout the paper. Continuous-time regulators will be denoted $K\left(s\right)\equiv K$, while their discrete counterparts, as a function of the sampling period will be used as $K_T\left(z\right)\equiv K_T$. The set of continuous-time systems will be denoted by $\mathcal{G}$, while the set of discrete-time systems will be denoted by ${\mathcal{G}}_D$. The standard growth functions and their notations according to [23]:

$$\begin{array}{cc}\hfill \mathcal{O}\left(g\right(n\left)\right)& =\left\{f\left(n\right)\exists c,n_0>0\mathrm{such}\mathrm{that}0\le f\left(n\right)\le cg\left(n\right),\forall n\ge n_0\right\};\hfill \\ \hfill \mathrm{\Omega }\left(g\right(n\left)\right)& =\left\{f\left(n\right)\exists c,n_0>0\mathrm{such}\mathrm{that}0\le cg\left(n\right)\le f\left(n\right),\forall n\ge n_0\right\};\hfill \\ \hfill \mathrm{\Theta }\left(g\right(n\left)\right)& =\left\{f\left(n\right)\exists c_1,c_2,n_0>0\mathrm{such}\mathrm{that}0\le c_{1}g\left(n\right)\le f\left(n\right)\le c_{2}g\left(n\right),\forall n\ge n_0\right\}.\hfill \end{array}$$

2. Mathematical Background and Proposed Extensions

2.1. Numerical Closed-Loop Control

Modern control systems make use of numerically implemented regulators to coordinate continuous-time processes. The classical configuration of a numerical regulator is presented in Figure 1. It uses sample and hold circuits as interfaces to the continuous-time adjacent components and, as such, it imposes the zero-order hold discretization method for the plant G. As such, define the plant discretization methodology as the mapping $\left(G\left(s\right),T\right)\mapsto G_T\left(z\right)\equiv G\left(z\right)\in {\mathcal{G}}_D$:

$$G_T\left(z\right)=\mathcal{Z}\left\{{\mathcal{L}}^{-1}\left\{G_{zoh,T}\left(s\right)·G\left(s\right)\right\}\right\}:\mathcal{G}\times {\mathbb{R}}_{+}\to {\mathcal{G}}_D.$$

(1)

with $\mathcal{G}$ denoting the continuous-time model set, ${\mathcal{G}}_D$ being the discrete-time model set, and with the digital-to-analog converter, i.e., zero-order hold, model:

$$G_{zoh,T}\left(s\right)=\frac{1-e^{-sT}}{s}.$$

(2)

For the purpose of this paper, we consider a two-degrees-of-freedom (2DOF) control structure as in Figure 2, which has an inner controller $K_{\mathrm{in}}\left(s\right)$ usually designed for disturbance rejection, and a feedforward controller $K_{\mathrm{ff}}\left(s\right)$ which provides the necessary compensation for the steady-state tracking behavior, along with the process model with disturbance, described by $G\left(s\right)$ and $G_d\left(s\right)$. As $K_{\mathrm{ff}}$ does not influence the signal path from the disturbance $d\left(t\right)$ to the output $y\left(t\right)$, the usual design workflow is to synthesize $K_{\mathrm{in}}$ and, then, to fine-tune the transient response from $r\left(t\right)$ to $y\left(t\right)$ through $K_{\mathrm{ff}}$. Both continuous-time controllers $K_{\mathrm{in}}\left(s\right)$ and $K_{\mathrm{ff}}\left(s\right)$ are assumed to have the following linear and time-invariant (LTI) state-space representations:

$$\left(K_x\left(s\right)\right):\left\{\begin{array}{cc}{\dot{\mathbf{x}}}_c\hfill & =A_{K,x}{\mathbf{x}}_c+B_{K,x}{\mathbf{u}}_c\hfill \\ {\mathbf{y}}_c\hfill & =C_{K,x}{\mathbf{x}}_c+D_{K,x}{\mathbf{u}}_c\hfill \end{array}\right.,x\in \{\mathrm{in},\mathrm{ff}\}.$$

(3)

Given that the continuous-time controllers $K_{\mathrm{in}}\left(s\right)$ and $K_{\mathrm{ff}}\left(s\right)$ must be numerically implemented on a microcontroller, the problem of selecting an appropriate sampling period $T\in (0,\infty )$ becomes a critical step. For a fixed sampling rate T, the discrete-time form of a controller $K_x\left(s\right)$ can be written as:

$$\left(K_x\left(z\right)\right):\left\{\begin{array}{cc}{\mathrm{x}}_{k+1}\hfill & ={\mathrm{\Phi }}_x{\mathrm{x}}_k+{\mathrm{\Gamma }}_x{\mathbf{u}}_k\hfill \\ {\mathbf{y}}_k\hfill & =C_{K,x}{\mathrm{x}}_k+D_{K,x}{\mathbf{u}}_k\hfill \end{array}\right.,x\in \{\mathrm{in},\mathrm{ff}\}.$$

(4)

To compute the transition matrix ${\mathrm{\Phi }}_x$ and the input matrix ${\mathrm{\Gamma }}_x$, a wide selection of discretization methods, denoted by $\mathbf{D}$, can be considered, such as zero-order hold, trapezoidal (Tustin), forward or backward Euler, frequency-response regression and so on. Additionally, in most situations, the output matrix remains the same as in the continuous-time representation, but there exist cases when it may be different than its continuous-time counterpart. As such, from a tuple of a continuous-time transfer matrix $K_x\left(s\right)=(A_{K,x},B_{K,x},C_{K,x},D_{K,x})\in \mathcal{G}$ and a sampling rate $T\in {\mathbb{R}}_{+}$ results an equivalent discrete-time transfer matrix $\left(K\left(s\right),T\right)\mapsto K_T\left(z\right)\equiv K\left(z\right)\in {\mathcal{G}}_D$, where:

$$K_T\left(z\right)=\mathbf{D}\left\{K\left(s\right),T\right\}:\mathcal{G}\times {\mathbb{R}}_{+}\to {\mathcal{G}}_D.$$

(5)

The open-loop process model has the input-output representation:

$$Y\left(s\right)=G_d\left(s\right)·D\left(s\right)+G\left(s\right)·U\left(s\right),$$

(6)

while the closed-loop system expression becomes:

$$Y\left(s\right)=\frac{G_d\left(s\right)}{1+G\left(s\right)K_{\mathrm{in}}\left(s\right)}·D\left(s\right)+\frac{G\left(s\right)(K_{\mathrm{in}}\left(s\right)+K_{\mathrm{ff}}\left(s\right))}{1+G\left(s\right)K_{\mathrm{in}}\left(s\right)}·R\left(s\right).$$

(7)

The discrete-time equivalent closed-loop model will be:

$$Y\left(z\right)=\frac{G_d\left(z\right)}{1+G\left(z\right)K_{\mathrm{in}}\left(z\right)}·D\left(z\right)+\frac{G\left(z\right)(K_{\mathrm{in}}\left(z\right)+K_{\mathrm{ff}}\left(z\right))}{1+G\left(z\right)K_{\mathrm{in}}\left(z\right)}·R\left(z\right)=H_{cl}^d\left(z\right)D\left(z\right)+H_{cl}^r\left(z\right)R\left(z\right),$$

(8)

which represents the starting point for various control design methodologies.

2.2. Sampling Rate Optimization

The current subsection presents an extension of the work from the paper [21].

2.2.1. Linear System Sampling Background

The well-known Shannon theorem [1] states that the necessary sampling rate to avoid information loss must be at most twice time smaller than the inverse of the continuous-time signal’s maximum frequency component:

$$T_s<\frac{T_{\min }}{2}\iff f_s>2·f_{\max },\mathrm{where}{f}_{\max }=\frac{1}{T_{\min }}.$$

(9)

Moreover, the theoretical upper bound of the frequency range after which aliasing phenomena occur is called Nyquist frequency and is given by:

$${\omega }_N=\frac{{\omega }_s}{2}=\pi ·f_s.$$

(10)

The above-mentioned version of Nyquist-Shannon theorem for signals can be extended to LTI systems by considering that the resulting discrete-time system must maintain all relevant dynamics of the continuous-time system. The dynamics of a system is given by its set of poles, one possible constraint for the sampling rate implying to be at most twice smaller than the inverse of the real part of the minimum value of poles. However, there are cases when the imaginary part of the poles is relevant, or even the values of the transmission zeros, for a good representation of the system’s frequency response.

For the remainder of the paper, the sampling rate will be denoted by $T_s\equiv T$. The analytical relationship between the continuous-time s-plane and the discrete-time z-plane is illustrated by the equation:

$$z=e^{s·T}.$$

(11)

Considering a point $s=\sigma +j\omega$, with $\sigma ,\omega \in \mathbb{R}$, the resulting corresponding complex number is $z=e^{\sigma T}\left(\cos \left(\omega T\right)+jsin\left(\omega T\right)\right)$. As such, the left half of the complex s-plane can be split into an infinite number of disjoint strips of height $\frac{2\pi }{T}$, due to the periodical nature of the $sin(·)$ and $cos(·)$ functions. The resulting primary strip contains points which correspond to $\omega \in \left[-{\omega }_N,{\omega }_N\right]$. The s-plane to z-plane mapping causes the following topologies of the primary strip of the s-plane:

• the imaginary axis in is mapped to the unit circle $\mathcal{C}\left(0,1\right)$;
• the upper and the lower edges are both mapped to the negative axis;
• the negative real axis is mapped to the positive axis inside the unit circle;
• the interior of the primary strip is mapped in the unit disk.

When sampling an LTI controller, several remarks can be considered then the behavior of the singularities is studied:

(i)

unstable zeros and poles tend decreasingly to the unit circle’s circumference, i.e., $\left|z\right|\searrow 1$ and $\left|p\right|\searrow 1$ as $T\searrow 0$; for poles, this aspect is relevant when sampling systems, as controllers generally do not employ unstable poles in their structure;

(ii)

stable zeros and poles tend increasingly to the unit circle’s circumference, i.e., $\left|z\right|\nearrow 1$ and $\left|p\right|\nearrow 1$ as $T\searrow 0$, requiring additional decimal or binary digits for an accurate representation.

(iii)

poles and zeros from the imaginary axis in the s-domain are maintained on the unit circle in the z-domain irrespective of the sampling period $T>0$. Additionally, integrator and derivative terms do not matter in deciding the practical sampling time.

(iv)

the quantization error increases as $T\searrow 0$ in the case of stable closed-loop systems.

Remark 1.

The quantization error mentioned in (iv) increases as $T\searrow 0$ as proved in [24], due to the guaranteed steady-state error being proportional to $\frac{1}{1-\rho \left(\mathrm{\Phi }\right)}$, while $\rho \left(\mathrm{\Phi }\right)$ denotes the spectral radius of the closed-loop state matrix: $\mathrm{\Phi }=\widehat{A}-\widehat{B}\widehat{D}\widehat{C}$, based on the series interconnection between the controller and the process model matrices.

According to the previously mentioned aspects, a sub-optimal value of the sampling rate must ensure a good trade-off between the fidelity of the representation of the continuous-time controllers and their implementability difficulty. As such, we next introduce a set of functionals to quantify these two aspects.

2.2.2. Proposed Functionals

To measure the similarity between a continuous-time LTI system $H\left(s\right)$ and a discrete-time system $H_T\left(z\right)$ over the frequency range $\mathrm{\Omega }$, the following functional can be defined as ${\mathcal{S}}_H^{\mathrm{\Omega }}:\mathcal{D}\to (0,\infty )$, with $\mathcal{D}=\mathcal{G}\times {\mathbb{R}}_{+}$:

$${\mathcal{S}}_H^{\mathrm{\Omega }}\left(T\right)={\int }_{\mathrm{\Omega }}{\left(\overline{\sigma }\left(H\left(j\omega \right)\right)-\overline{\sigma }\left(H_T\left(e^{j\omega T}\right)\right)\right)}^2\mathrm{d}\omega ,$$

(12)

where $\overline{\sigma }(·)$ is the maximum singular value. Moreover, as mentioned before, the available frequency domain is $\mathrm{\Omega }=(0,{\omega }_N)$, ${\omega }_N$ being the corresponding Nyquist-Shannon frequency to the sampling rate T. Therefore, the similarity functional $S_H:\mathcal{D}\to (0,\infty )$ becomes:

$${\mathcal{S}}_H\left(T\right)={\int }_{0^{+}}^{{\omega }_N^{-}}{\left(\overline{\sigma }\left(H\left(j\omega \right)\right)-\overline{\sigma }\left(H_T\left(e^{j\omega T}\right)\right)\right)}^2\mathrm{d}\omega .$$

(13)

However, to compute the similarity integral term ${\mathcal{S}}_H$, a discrete set from the frequency range $(0,{\omega }_N)\to \overline{\mathrm{\Omega }}=\left\{\underline{\omega }={\omega }_1<{\omega }_2<\cdots <{\omega }_N-{\omega }_{\epsilon }=\overline{\omega }\right\}$ can be considered, the performance index being approximated as follows:

$${\mathcal{S}}_H^{\overline{\mathrm{\Omega }}}\left(T\right)\simeq \sum _{\omega \in \overline{\mathrm{\Omega }}}{\left(\overline{\sigma }\left(H\left(j\omega \right)\right)-\overline{\sigma }\left(H_T\left(e^{j\omega T}\right)\right)\right)}^2\Delta \omega ,$$

(14)

where ${\omega }_{\epsilon }>0$ is a predefined threshold used to avoid the prewarping phenomenon which ensues in the magnitude responses when $\omega \to {\omega }_N$.

Remark 2.

There are several alternatives in defining the similarity between two LTI systems, with particular interest for closed-loop connections, such as the normalized dissimilarity function metric, described in [25], the standard ν-gap metric, with its limitations exposed in [26], or the improved Vinnicombe ν-gap metric [27], described as:

$$\delta (G_1,G_2)=max\left(\overline{\delta }\left(G_1,G_2\right),\overline{\delta }\left(G_2,G_1\right)\right),$$

(15)

with:

$$\overline{\delta }(G_1,G_2)=\underset{Q\in H_{\infty }}{inf}{∥\left(\begin{array}{c}{M}_1\\ N_1\end{array}\right)-\left(\begin{array}{c}{M}_2\\ N_2\end{array}\right)Q∥}_{\infty },$$

(16)

where for the reference systems $G_1$ and $G_2$, their right normalized coprime factorizations will be considered: $G_1=N_1{M}_1^{-1}$ and $G_2=N_2{M}_2^{-1}$, which in context of the proposed sampling rate optimization problem, the differences should be computed with the subsystems applied in $G_1\equiv K\left(s\right)\mapsto j\omega$ and $G_2\equiv K_T\left(z\right)\mapsto e^{j\omega T}$, respectively.

To define the fixed-point implementability functional ${\mathcal{I}}_H\left(T\right):\mathcal{D}\to (0,\infty )$ which measures the quantization implementation difficulty of the LTI system $H_T\left(z\right)$, the following expression can be considered:

$${\mathcal{I}}_H\left(T\right)=\frac{1}{min\left\{\left|\left|\lambda \right|-1\right|,\lambda \in \mathcal{P}\left(H_T\right)\cup \mathcal{Z}\left(H_T\right)\right\}},$$

(17)

with $\mathcal{P}$ and $\mathcal{Z}$ denoting the pole and zero sets of $H_T$, respectively, excluding singularities from the set $\mathcal{C}\left(0,1\right)$.

Alongside the previous implementability functional, an execution time cost functional ${\mathcal{T}}_H\left(T\right):\mathcal{D}\to (0,\infty )$ becomes necessary to limit the decrease of $T\searrow 0$ which would ideally lead the similarity functional costs to zero:

$${\mathcal{T}}_H\left(T\right)=\frac{1}{T}.$$

(18)

A concluding functional with global effect, ${\mathcal{J}}_{\mathrm{stab}\left(H\right)}:\mathcal{D}\to \{0,\infty \}$, will define the feasibility domain of the optimization problem, because it accepts or rejects a specific sampling period value. It induces an infinitely valued constant when the numeric system $H_T$ becomes unstable, and otherwise, defaults to no extra penalization:

$${\mathcal{J}}_{\mathrm{stab}\left(H\right)}\left(T\right)=\left\{\begin{array}{c}+\infty ,\mathrm{if}{H}_T\mathrm{has}\mathrm{unstable}\mathrm{poles};\hfill \\ 0,\mathrm{otherwise}.\hfill \end{array}\right.$$

(19)

For the numerical implementation of the stability functional, a sufficiently large value ${\alpha }_{\infty }$ can be considered to mark the feasibility subdomain of $\mathcal{D}$ where the system H is stable, leading to an approximation of the original performance index:

$${\mathcal{J}}_{\mathrm{stab}\left(H\right)}^{{\alpha }_{\infty }}\left(T\right)=\left\{\begin{array}{c}+{\alpha }_{\infty },\mathrm{if}{H}_T\mathrm{has}\mathrm{unstable}\mathrm{poles};\hfill \\ 0,\mathrm{otherwise}.\hfill \end{array}\right.$$

(20)

2.2.3. Optimization Problem

For the two-degrees-of-freedom control structure proposed in Figure 2 the following functionals will be considered to formulate the optimization problem to determine the sampling rate:

• two terms ${\mathcal{S}}_{K_{\mathrm{in}}}\left(T\right)$ and ${\mathcal{S}}_{K_{\mathrm{ff}}}\left(T\right)$ representing the similarity between the continuous-time and the discrete-time representations of the controllers $K_{\mathrm{in}}\left(s\right)$ and $K_{\mathrm{ff}\left(s\right)}$ over the frequency range $(0,{\omega }_N)$;
• two terms ${\mathcal{S}}_{H_{cl}^r}\left(T\right)$ and ${\mathcal{S}}_{H_{cl}^d}\left(T\right)$ representing the similarity between the continuous-time and the discrete-time representations of the resulting closed-loop systems $H_{cl}^r\left(s\right)$ and $H_{cl}^d\left(s\right)$ over the same frequency range;
• the quantization implementation difficulty ${\mathcal{I}}_{K_{\mathrm{in}}}\left(T\right)$ and ${\mathcal{I}}_{K_{\mathrm{ff}}}\left(T\right)$ of each controller along with the execution time functional $\mathcal{T}\left(T\right)$;
• the stability functional ${\mathcal{J}}_{\mathrm{stab}\left(H_{cl}^r\right)}\left(T\right)$ of the resulting numerical closed-loop system $H_{cl,T}^r\left(z\right)$.

The first set of three terms are used to measure the fidelity between the components of the continuous-time designed system and the resulting components in the discrete-time domain, considering the transient and steady-state performance altering. The next set of three terms manages to encompass the implementation difficulty of the resulting controller, according to the already-mentioned remarks. The last term is nothing but a correction used to define the feasibility subdomain of $\mathcal{D}$. As such, a trade-off between the first two sets of performance indices must be established by taking into account the last feasibility term, resulting a non-convex optimization problem. However, in order to increase the flexibility of this optimization problem, a set of weights $c_{1-6}$ can be considered as follows: the weights $c_1$, $c_2$, $c_3$, and $c_4$ are for the fidelity measurement indices, while weights $c_5$, $c_6$, and $c_7$ are for the implementability indices, the final functional $\mathcal{J}:\mathcal{D}\to {\overline{\mathbb{R}}}_{+}$ being:

$$\mathcal{J}\left(T\right)=c_1{\mathcal{S}}_{K_{\mathrm{in}}}\left(T\right)+c_2{\mathcal{S}}_{K_{\mathrm{ff}}}\left(T\right)+c_3{\mathcal{S}}_{H_{cl}^r}\left(T\right)+c_4{\mathcal{S}}_{H_{cl}^d}\left(T\right)c_5{\mathcal{I}}_{K_{\mathrm{in}}}\left(T\right)+c_6{\mathcal{I}}_{K_{\mathrm{ff}}}\left(T\right)+c_7\mathcal{T}\left(T\right)+{\mathcal{J}}_{\mathrm{stab}\left(H_0\right)}\left(T\right).$$

(21)

Moreover, considering the numerical approximations (14) and (20), the following numerically implementable non-convex optimization problem occurs:

Problem 1.

For an LTI plant model $G\left(s\right)$ included in a two-degrees-of-freedom control structure with the continuous-time inner controller $K_{\mathrm{in}}\left(s\right)$ and with the continuous-time feedforward controller $K_{\mathit{ff}}\left(s\right)$ as in Figure 2, define an ordered set of pulsations, preferably in logarithmic scale:

$$\overline{\mathrm{\Omega }}=\left(\underline{\omega }={\omega }_1<{\omega }_2<\dots <{\omega }_N-{\omega }_{\epsilon }=\overline{\omega }\right).$$

(22)

Then, the functional ${\mathcal{J}}^{\overline{\mathrm{\Omega }}}:\mathcal{D}\to {\mathbb{R}}_{+}$:

$${\mathcal{J}}^{\overline{\mathrm{\Omega }}}\left(T\right)=c_1{\mathcal{S}}_{K_{\mathrm{in}}}^{\overline{\mathrm{\Omega }}}\left(T\right)+c_2{\mathcal{S}}_{K_{\mathrm{ff}}}^{\overline{\mathrm{\Omega }}}\left(T\right)c_3{\mathcal{S}}_{H_{cl}^r}^{\overline{\mathrm{\Omega }}}\left(T\right)+c_4{\mathcal{S}}_{H_{cl}^d}^{\overline{\mathrm{\Omega }}}\left(T\right)+c_5{\mathcal{I}}_{K_{\mathrm{in}}}\left(T\right)+c_6{\mathcal{I}}_{K_{\mathrm{ff}}}\left(T\right)+c_7\mathcal{T}\left(T\right)+{\mathcal{J}}_{\mathrm{stab}\left(H_0\right)}^{{\alpha }_{\infty }}\left(T\right),$$

(23)

with weighting terms $c_{1-7}$, leads to the following quasi-optimal sampling time optimization problem:

$$\underset{T\in \mathcal{D}}{min}{\mathcal{J}}^{\overline{\mathrm{\Omega }}}\left(T\right).$$

(24)

Remark 3.

The resulting optimization problem (24) is non-convex by nature. As such, the line-search procedure used to solve this problem must be able to explore the whole feasible domain. One possible solution is a metaheurisitc approach, such as particle swarm optimization (PSO) [28], which will be used for the purpose of this paper.

2.3. Worst-Case Execution Time Analysis

The current subsection presents an extension of the work from the paper [10].

2.3.1. Execution Time Model

This subsection proposes the study of the implementation details for the case of linear control structures, where controllers will be modeled as infinite-impulse response (IIR) and finite-impulse response (FIR) filters. Starting from [10], a set of implementation aspects which are treated for state-space realizations will be considered in a unified manner for the case of transfer matrices in this paper. Additionally, this section also analyzes the 2DOF extension. The duration for each operation involved in the control structure can significantly impact the regulator implementability.

The necessary mathematical operations to fully implement an LTI-based control law must be formally defined. Besides the LTI-control law, the classical saturation and anti-windup nonlinearities will be considered, which are usually related to said LTI laws. The unary and binary mathematical operators defined in

Table 1. Formal operators necessary for the implementation of LTI-based control laws.

#	Operator	Domain	Definition	Observations
1	Addition	$\mathrm{𝒶}:{\mathbb{R}}^2\to \mathbb{R}$	$\mathrm{𝒶}\left(x_1,x_2\right)=x_1+x_2$	bilinear
2	Multiplication	$\mathrm{𝓂}:{\mathbb{R}}^2\to \mathbb{R}$	$\mathrm{𝓂}\left(g,x\right)=g·x$	homogenous
3	Saturation	${\mathrm{𝓈}}_{\underline{x},\overline{x}}:\mathbb{R}\to \mathbb{R}$	${\mathrm{𝓈}}_{\underline{x},\overline{x}}\left(x\right)=\left\{\begin{array}{cc}\underline{x},\hfill & \mathrm{if}x<\underline{x};\hfill \\ x,\hfill & \mathrm{if}\underline{x}\le x\le \overline{x};\hfill \\ \overline{x},\hfill & \mathrm{if}\overline{x}<x,\hfill \end{array}\right.$	nonlinear
4	Anti-windup	${\mathrm{𝓌}}_{\underline{x},\overline{x}}:\mathbb{R}\to \mathbb{R}$	various [12,29,30]	for integrators
5	Load/Store	$\mathrm{𝓁}\left(x\right):\mathbb{R}\to \mathbb{R}$	$\mathrm{𝓁}\left(x\right)=x$	bilinear
6	Null	$\mathrm{𝓃}\left(x\right):\mathbb{R}\to \mathbb{R}$	$\mathrm{𝓃}\left(x\right)=0$	for delays

and gathered in the operations alphabet $\mathcal{O}=\left\{\mathrm{𝓃},\mathrm{𝒶},\mathrm{𝓂},\mathrm{𝓈},\mathrm{𝓌},\mathrm{𝓁}\right\}$, are considered with real operands and must be accounted for into a microprocessor-based environment.

As such, the process of computing a command signal $\mathbf{y}\left[k\right]$ as in Figure 1 implies a finite and formal computational finite sequence $\mathcal{S}𝒸\in {\mathcal{O}}^N$, where all terms are mathematical operations as in Table 1, i.e., $\mathcal{S}𝒸\left[i\right]=p_i\in \mathcal{O}$, $i=\overline{1,N}$, where N depends on the structure of the controller $K\left(z\right)$. Moreover, in order to additionally specify a set of practical hardware specifications and constraints $\mathcal{H}$ and to uniformly describe the problem, the finite sequence can be now extended to a full sequence ${\mathcal{S}𝒸}^{K,\mathcal{H}}\in {\mathcal{O}}^{\infty }$:

$${\mathcal{S}𝒸}^{K,\mathcal{H}}=\left(p_1,p_2,\cdots ,p_N,\mathcal{n},\mathcal{n},\cdots \right),p_i\in \mathcal{O}.$$

(25)

Starting from an array ${\mathcal{S}𝒸}^{K,\mathcal{H}}$ of operations as (25), we follow with a general-purpose instruction set model. Assume a Random-Access Machine (RAM) computational model as in [23], with deterministic operations. RAM machines have practical counterparts, materialized through RISC machines. Reconfigurable RISC machines specialized on certain problems have been proposed in [31]. Additionally, there is the approach of multiply and accumulate (MAC) instructions supported in digital signal processors (DSP) [32]. Depending on the supported computer architectures of the RCP framework, relevant are also Single Instruction stream/Multiple instruction Pipelining (SIMP) [33] constructions with respect to single-processor architectures, or Single Instruction/Multiple Data (SIMD) features, which allow the practical parallelization of addition and multiplication operations for several sets of operands.

Definition 1.

The hardware constraints and specification set $\mathcal{H}$ encompasses metadata which imply extra operations or different approaches to the standard operations performed on the controller signals and implementation-specific information, with various outcomes on the total execution time, frequently found in practice being:

• reading reference signals $r\left[k\right]$ and plant measurements $y\left[k\right]$, all input signal reading steps may imply preprocessing constraints in terms of sensor delays, impulse counters or data type conversions—this equates to adding $p_i\in \left\{\mathrm{𝓃},\mathrm{𝓁}\right\}$ steps;
• scaling operations for the input and output signals imposed by the operating point used for plant linearization: $\Delta u\left[k\right]=u\left[k\right]-u_0$ and $y\left[k\right]=y_0+\Delta y\left[k\right]$, which equates to augmenting the sequence set with $p_i\in \left\{\mathrm{𝒶},\mathrm{𝓈},\mathrm{𝓁}\right\}$ items;
• input and output signals scaling operations, i.e., $u_s\left[k\right]=a_u·u\left[k\right]+b_u$ and $y_s\left[k\right]=a_y·y\left[k\right]+b_y$, useful especially for sensor/adapter signals, and extend the operations with $p_i\in \left\{\mathrm{𝒶},\mathrm{𝓂},\mathrm{𝓈},\mathrm{𝓁}\right\}$;
• starting from the variable base word length L of the microcontroller arithmetic registers which allows operations to be executed in a single clock tick, each variable’s type and size should be adequately adapted for $2\times L$, $4\times L$ etc., which complicates the adding and multiplication routines with additional $p_i\in \left\{\mathrm{𝒶},\mathrm{𝓂},\mathrm{𝓁}\right\}$ steps;
• controller gain-scheduling verifications and updates based on the value of the input signal $u\left[k\right]$, leading to extra $p_i\in \left\{\mathrm{𝓁}\right\}$;
• underflow and overflow checks for involved signals, implying saturations $p_i\in \left\{\mathrm{𝓈}\right\}$;
• availability of Direct Memory Access (DMA) modules, MAC instructions, circular buffers in opposition to linear buffering, or output bypassing, which has the advantage of discarding $p_i\in \left\{\mathrm{𝓁}\right\}$ operations.

Such specifications will be quantified by scaling factors ${\gamma }_j$ to the base duration of the RISC machine model operations. Depending on the significance of the alternative instruction, ${\gamma }_j$ could be less than, equal to, or greater than one, respectively.

Therefore, a set of instructions ${\mathcal{S}𝒸}^{K,\mathcal{H}}$ can be manually designed or automatically deduced. For manual modeling, the control engineer needs to find this set for each control law, while for automatic mode, an RCP tool can already deduce this set. Such an RCP tool generally deals with the code generation for different environments. Now, the sequence of operations generator procedures can be seen as a functional:

$$\mathrm{\Psi }:{\mathcal{G}}_D\times \mathcal{H}\to {\mathcal{O}}^{\infty },\mathrm{\Psi }\left(K,\mathcal{H}\right)={\mathcal{S}𝒸}^{K,\mathcal{H}}.$$

(26)

To implement the mathematical operations $p_i\in \mathcal{O}$ from Table 1, the atomic assembly instructions $a_i\in \mathcal{A}=\left\{\mathtt{NOP},\mathtt{MF},\mathtt{MS},\mathtt{ADD},\mathtt{MUL},\mathtt{SH},\mathtt{JMP},\mathtt{CMP}\right\}$ from

Table 2. Base assembly operations $a_i\in \mathcal{A}$ for a Random-Access Machine model used in an LTI-based control system context.

#	Operator	Abbreviation	Part of Operation	Observations
1	No operation	NOP	$\mathrm{𝓃}$	–
2	Memory fetch	MF	$\mathrm{𝓁}$	SIMP
3	Memory store	MS	$\mathrm{𝓁}$	SIMP
4	Add	ADD	$\mathrm{𝒶}$	SIMP, MAC, SIMD
5	Multiply	MUL	$\mathrm{𝓂}$	SIMP, MAC, SIMD
6	Binary shift	SH	$\mathrm{𝓂}$	SIMP, SIMD
7	Jump	JMP	$\mathrm{𝓈},\mathrm{𝓌},\mathrm{𝓁}$	–
8	Compare	CMP	$\mathrm{𝓈},\mathrm{𝓌},\mathrm{𝓁}$	SIMP

will represent the starting point. It covers the basic arithmetical operations required in linear systems, with the additivity and homogeneity properties, conditional jumps for saturations and the anti-windup of integrator terms, forced and imposed delays through data acquisition hardware and access to memory devices. Each arbitrary atomic assembly operation will be denoted by a, as part of the formal set $\mathcal{A}$ all the equivalent assembly instructions supported by $\mathcal{H}$ for the implementation of $\mathcal{S}𝒸$, with the RISC machine assumption that each instruction takes a fixed clock tick value $T_{clk}$ > 0. The number of ways in obtaining the resulting assembly instructions is not unique and it further depends on the structure of $\mathcal{H}$. A straightforward example to illustrate this phenomenon is when the regulator coefficients are stored contiguously in the memory in comparison to arbitrary and uncorrelated memory registers. The execution time implications are obtained through different pointer operations. To conclude, a new mathematical operator analogous to (26), tasked with the generation of a computer-equivalent set of instructions given by the functional $\mathcal{S}𝓅\in {\mathcal{A}}^{\infty }$ is:

$$\mathrm{\Xi }:{\mathcal{G}}_D\times \mathcal{H}\to {\mathcal{A}}^{\infty },\mathrm{\Xi }\left(K,\mathcal{H}\right)={\mathcal{S}𝓅}^{K,\mathcal{H}},$$

(27)

which results in an infinite sequence of atomic software instructions, but with a finite number of them being different to NOP, located at the start of the sequence, which implement the linear controller formula:

$${\mathcal{S}𝓅}^{K,\mathcal{H}}=\left(a_1,\cdots ,a_M,\mathtt{NOP},\mathtt{NOP},\cdots \right),a_i\in \mathcal{A}.$$

(28)

The difference between the sets $\mathcal{S}𝒸$ and $\mathcal{S}𝓅$ is that $\mathcal{S}𝒸$ contains abstract mathematical operations $p_i\in \mathcal{O}$, and each such operation $p_i$ will be practically implemented using equivalent $a_{i,j}\in \mathcal{A}$ steps, $j=\overline{1,N_i}$, as in the mapping:

$$p_i\mapsto \left(a_{i,1},a_{i,2},\cdots ,a_{i,N_i}\right),\forall i=\overline{1,N},$$

(29)

the number of atomic operations different from NOP being $M=N_1+N_2+\cdots +N_N$.

The importance of the previously defined sequences and functionals, i.e., ${\mathcal{S}𝒸}^{K,\mathcal{H}}$, $\mathsf{\Psi }$, ${\mathcal{S}𝓅}^{K,\mathcal{H}}$ and $\mathsf{\Xi }$, respectively, and the implications of estimating the number of assembly operations in a tight manner was insisted upon in the base paper [10]. Figure 3 gathers them and illustrates their connections in an RCP context. The mapping from $\left(K,\mathcal{H}\right)$ to the set ${\mathcal{S}𝒸}^{K,\mathcal{H}}$ is usually performed for Model-in-the-Loop (MiL) simulations through an application $\mathsf{\Psi }$, while the mapping $\left(K,\mathcal{H}\right)\mapsto {\mathcal{S}𝓅}^{K,\mathcal{H}}$ is made through Software-in-the-Loop (SiL) testing using an application $\mathsf{\Xi }$. The master RCP program, with access to both ${\mathcal{S}𝒸}^{K,\mathcal{H}}$ and ${\mathcal{S}𝓅}^{K,\mathcal{H}}$, can subsequently perform a worst-case execution time analysis on the implementation of the digital regulator $K\left(z\right)$ in the production hardware context $\mathcal{H}$.

Figure 4 presents the sequence diagram with the timing constraints of the discrete-time controller K with respect to other software threads from the microcontroller, with an illustration of the WCET of the controller interrupt service routine (ISR) thread. The main result of the section is gathered in the following theorem.

Theorem 1.

Given a numerical control law given by $K\in {\mathcal{G}}_D$, along with a microcontroller specification set $\mathcal{H}$ and a code-generation procedure given by a pair of MiL and SiL $\left(\mathrm{\Psi },\mathrm{\Xi }\right)$, the worst-case execution time estimation can be computed by the following formula:

$${\displaystyle \mathit{WCET}}\left(\mathrm{\Psi },\mathrm{\Xi }\right)=\left(\left|{\mathcal{S}𝓅}^{K,\mathcal{H}}\right|+\mathcal{O}\left(1\right)\right)\times T_{clk},$$

(30)

where $\left|{\mathcal{S}𝓅}^{K,\mathcal{H}}\right|$ represents the number of atomic assembly operations $a_i\in \mathcal{A}$ as in (28), and $\mathcal{O}\left(1\right)$ accounts the context switching operations for the other software threads. Additionally, the exact bounds from $\mathcal{O}\left(1\right)$ depend on all other software entities running on the same microprocessor and are not correlated with the input dimension m or the output dimension p of the numerical controller.

Proof. 

According to the schematic representation from Figure 3, starting from the control law $K\in {\mathcal{G}}_D$ and the specification set $\mathcal{H}$, the set ${\mathcal{S}𝒸}^{K,\mathcal{H}}$ can be obtained via the MiL operator $\mathrm{\Psi }$, resulting in $p_1,p_2,\cdots ,p_N\in \mathcal{O}$. However, in order to measure each $p_i$, the set of atomic operations must be attached via the SiL operator $\mathrm{\Xi }$: $p_i\mapsto \left(a_{i,1},a_{i,2},\cdots ,a_{i,N_i}\right)$, each such atomic operation requiring exactly $T_{clk}$, leading to:

$$T_1\equiv t\left(K\right)=\sum _{i=1}^N\sum _{l=1}^{N_i}t\left(a_{i,l}\right)=T_{clk}·\sum _{i=1}^N\sum _{l=1}^{N_i}1=T_{clk}\times \left|{\mathcal{S}𝓅}^{K,\mathcal{H}}\right|,$$

(31)

where $t(·)$ is the time necessary to execute an operation or a set of operations.

Additionally, a second term $T_2={\sum }_{\mathrm{ISR}}\lambda$ accumulates ISR switching and stack-handling operations handled by the scheduler, bounded by a processor-dependent constant $\lambda >0$, which can be modeled as $\mathcal{O}\left(1\right)\times T_{clk}$. As such, the worst-case execution time can now be written as:

$$\mathtt{WCET}\left(\mathrm{\Psi },\mathrm{\Xi }\right)=T_1+T_2=\left(\left|{\mathcal{S}𝓅}^{K,\mathcal{H}}\right|+\mathcal{O}\left(1\right)\right)\times T_{clk},$$

(32)

which concludes the proof.    □

Observation 1.

All possible delays caused by context switching to preemptive ISRs belonging to measurement data processing, with the cost of $\mathcal{O}\left(m\right)$, are included in the input processing step and do not remain unaccounted for in the execution time model of Theorem (30).

Two additional performance qualifiers can be employed to globally assess the controller ISR implementation impact on the scheduling algorithm of the processor.

Definition 2.

The processor usage level qualifier relative to a fixed sampling period $T>0$ of a discrete-time regulator $K\left(z\right)\in {\mathcal{G}}_D$, described in a relative manner, is defined by:

$$\mathcal{U}\left(\mathrm{\Psi },\mathrm{\Xi },T\right)=\frac{{\displaystyle \mathit{WCET}}\left(\mathrm{\Psi },\mathrm{\Xi }\right)}{T}\times 100[\%].$$

(33)

Definition 3.

The processor idle time qualifier with respect to a fixed sampling period $T>0$ of a discrete-time regulator $K\left(z\right)\in {\mathcal{G}}_D$, described in absolute units, is defined by:

$$\mathcal{I}\left(\mathrm{\Psi },\mathrm{\Xi },T\right)=max\left\{0,T-{\displaystyle \mathit{WCET}}\left(\mathrm{\Psi },\mathrm{\Xi }\right)\right\}\left[s\right].$$

(34)

2.3.2. Modeling Duration of Finite and Infinite-Impulse Response Topologies

Denote by $\mathbf{H}\in {\mathcal{G}}_D^{p\times m}$ a MIMO regulator with m inputs and p outputs, thus fully described by the expressions of $m\times p$ transfer functions $H_{ij}\in {\mathcal{G}}_D$:

$$\left(\mathbf{H}\right):\left(\begin{array}{c}{Y}_1\left(z\right)\\ Y_2\left(z\right)\\ ⋮\\ Y_p\left(z\right)\end{array}\right)=\left(\begin{array}{cccc}{H}_{11}\left(z\right)& H_{12}\left(z\right)& \cdots & H_{1m}\left(z\right)\\ H_{21}\left(z\right)& H_{22}\left(z\right)& \cdots & H_{2m}\left(z\right)\\ ⋮& \ddots & \cdots & ⋮\\ H_{p1}\left(z\right)& H_{p2}\left(z\right)& \cdots & H_{pm}\left(z\right)\end{array}\right)·\left(\begin{array}{c}{U}_1\left(z\right)\\ U_2\left(z\right)\\ ⋮\\ U_m\left(z\right)\end{array}\right).$$

(35)

Each element of the transfer matrix $\mathbf{H}$, i.e., $H_{ij}$, can be modeled as an infinite-impulse response (IIR) filter or as a finite-impulse response (FIR) filter. The case where $\mathbf{H}$ is modeled as a state-space representation is treated in the base conference paper [10], namely in Algorithm 1 and Table I, respectively.

For an arbitrary discrete-time IIR transfer function $H\left(z\right)$ of order n, define $\mathrm{\Omega }\left(H\right)$ as the pair:

$$\mathrm{\Omega }\left(H\right)=\left(n_2,n_1\right)=\left\{\begin{array}{c}\left(⌊\frac{n}{2}⌋,1\right),\mathrm{if}n\mathrm{is}\mathrm{odd};\hfill \\ \left(\frac{n}{2},0\right),\mathrm{if}n\mathrm{is}\mathrm{even}.\hfill \end{array}\right.$$

(36)

Using this notation, the transfer function $H\left(z\right)$ can be written using a series of second-order sections and an additional first-order component, if necessary, as:

$$H_{IIR}\left(z\right)={\left(\frac{b_{0,1}z+b_{0,0}}{z+a_{0,0}}\right)}^{n_1}·\prod _{i=1}^{n_2}\left(\frac{b_{i,2}{z}^2+b_{i,1}z+b_{i,0}}{z^2+a_{i,1}z+a_{i,0}}\right)={\left((H_1\left(z\right)\right)}^{n_1}·\prod _{i=1}^{n_2}{H}_{2,i}\left(z\right),$$

(37)

with the terms from each triplet $\left(b_{i,2},b_{i,1},b_{i,0}\right)$ not all null, $H_1$ known as a first-order section, and each $H_{2,i}$, with $i=\overline{1,n_2}$ being denoted in the literature as a second-order section (SOS).

Remark 4.

Second-order sections are usually described with an additional gain term multiplied separately to the output of the transfer function as:

$$H_2\left(z\right)=g·\frac{b_2{z}^2+b_{1}z+b_0}{z^2+a_{1}z+a_0},g\ne 0.$$

(38)

To not further complicate the notations, we will not explicitly write this gain term for every second-order section, but it will be implicitly considered in the implementation and execution time analysis.

There are multiple approaches of implementing digital biquadratic filters, a good example being the description from the monograph [34].

Table 3. Digital biquadratic topology implementations; all difference equations implement the same input-output second-order transfer function, but differ through the configurations of the state signals.

IIR SOS Topology	Difference Equation	$min\left|{\mathit{𝒮}\mathit{𝒸}}^{\mathit{K},\mathcal{H}}\right|$
Direct Form I (DFI)	$y\left[k\right]=g·\left({\sum }_0^2{b}_{i}u[k-i]-{\sum }_1^2{a}_{i}y[k-i]]\right).$	$4\mathrm{𝒶}$, $6\mathrm{𝓂}$, $12\mathrm{𝓁}$
Direct Form II (DFII)	$\left\{\begin{array}{c}y\left[k\right]=g·{\sum }_0^2\left(b_i·x[k-i]\right);\hfill \\ x\left[k\right]=u\left[k\right]-{\sum }_1^2\left(a_i·x[k-i]\right).\hfill \end{array}\right.$	$4\mathrm{𝒶}$, $6\mathrm{𝓂}$, $14\mathrm{𝓁}$
$\begin{array}{c}\mathrm{Transposed}\mathrm{Direct}\\ \mathrm{Form}\mathrm{I}\left(\mathrm{TDFI}\right)\end{array}$	$\left\{\begin{array}{c}y\left[k\right]=g·{\sum }_0^2\left(b_i·x\right)[k-i];\hfill \\ x\left[k\right]=u\left[k\right]-{\sum }_1^2\left(a_i·x\right)[k-i].\hfill \end{array}\right.$	$4\mathrm{𝒶}$, $6\mathrm{𝓂}$, $14\mathrm{𝓁}$
$\begin{array}{c}\mathrm{Transposed}\mathrm{Direct}\\ \mathrm{Form}\mathrm{II}\left(\mathrm{TDFII}\right)\end{array}$	$\left\{\begin{array}{c}y\left[k\right]=g·\left(b_{0}u\left[k\right]+x_1[k-1]\right);\hfill \\ x_1\left[k\right]=b_{1}u\left[k\right]-a_{1}y\left[k\right]+x_2[k-1];\hfill \\ x_2\left[k\right]=b_{2}u\left[k\right]-a_{2}y\left[k\right].\hfill \end{array}\right.$	$4\mathrm{𝒶}$, $6\mathrm{𝓂}$, $16\mathrm{𝓁}$

exposes the four usual topologies, which principally implement the same input-output transfer function, but with important differences regarding numerical stability when selecting fixed-point or floating-point implementations. These configurations are referred to as the canonical forms: Direct Form I, Direct Form II, Transposed Form I, Transposed Form II, which differ in the numerical properties of their implementations and in the number of necessary delay elements. As observed in the third column of the table, all biquad topologies are based on four additions, five multiplications, and a different number of load/store operations, depending on the definition of the internal state variables. Such implementation details are relevant when studying particularized structures, such as in the situations treated in [35,36].

The ISR model for IIR controller structures is based on the pseudocode exposed in Algorithm 1. hlBased on the specifications for $\mathcal{H}$, each line of the pseudocode will have a set of mandatory mathematical operations, along with optional operations, which will be accounted for in the execution time analysis model through different constant weights.

Algorithm 1 Infinite-impulse response (IIR) filter interrupt service routine (ISR)

Input:  $H_{ij}\left(z\right)$ as in (35), topology as in Table 3 Output:  Execution time profiler for routine iirFiltIsr 1: Construct software structures $H_1$, $H_{2,1}$–$H_{2,n_2}$ according to (37) and Table 3 2: Initialize delays involved in second-order sections to zero 3: procedureiirFiltIsr($\mathrm{\Omega }\left(H\right)$,$H_2\left[\right]$$\left[,H_1\right]$) 4:     Read and scale $u\left[k\right]$ from input device 5:     $\tilde{u}$←$u\left[k\right]$ 6:     for $i\leftarrow 1$ to $n_2$ do 7:         $\tilde{u}$←$H_2\left[i\right].\mathtt{call}\left(\tilde{u}\right)$  ▷ The output from $H_2[i-1]$ becomes input to $H_2\left[i\right]$ 8:     end for 9:     $\left[\tilde{u}\leftarrow H_1.\mathtt{call}\left(\tilde{u}\right)\right]$ 10:     $y\left[k\right]=\tilde{u}$ 11:     Scale $y\left[k\right]$ and write to output device 12: end procedure 13: functioncall(H, topology) ▷ First or second-order section call method for IIR filter 14:     Read input $u\left[k\right]$ 15:     Shift input delays $u\left[k\right]$, $u[k-1]$$,\left[u[k-2]\right]$ 16:     Compute $y\left[k\right]$ according to input topology as in Table 3 17:     Shift output delays $y\left[k\right]$, $y[k-1]$$,\left[y[k-2]\right]$ 18:     Scale by second-order section gain g 19:     return $y\left[k\right]$ 20: end function

In ref. [10], an analysis has been performed on the simplified case of a nth order IIR SISO transfer function in a series connection, described as:

$$H_{IIR}^s\left(z\right)=\frac{Y\left(z\right)}{U\left(z\right)}=\frac{b_m{z}^m+b_{m-1}{z}^{m-1}+\cdots +b_0}{z^n+a_{n-1}{z}^{n-1}+a_{n-2}{z}^{n-2}+\cdots +a_0},$$

(39)

which, by design, it can implicitly include a delay $z^{-n_d}$ by forcing the first $n_d$ coefficients $b_i$ to zero. Such a transfer function has its corresponding difference equation as:

$$y\left[k\right]=-\sum _{i=0}^{n-1}{a}_i·y[k-i]+\sum _{j=0}^{m-1}{b}_j·u[k-j].$$

(40)

A further particularization on the structure of $H_{ij}\left(z\right)$ is to consider the expression of a FIR filter topology, which, by design, discards the previous output delays. The present command signal $\mathbf{y}\left[k\right]$ will, as such, depend only on an array of delays, i.e., delay tap of inputs $u[k-i]$, modeled as:

$$H_{FIR}\left(z\right)=g·z^{-n_d}·\sum _{i=0}^m{b}_i·z^{-i}\equiv g·\sum _{k=0}^{N-1}{h}_k·z^{-k},g\ne 0.$$

(41)

Four typical canonical forms are distinguished for FIR-type filters, namely the Direct Form (DF), Direct Form Transposed (DFT), Symmetric and Antisymmetric [34], respectively, with their definitions exposed in

Table 4. Digital FIR filter topology implementations; all difference equations implement the same input-output Nth order transfer function, but vary through the configuration of the accumulator.

FIR Topology	Difference Equation	$min\left|{\mathcal{S}\mathit{𝒸}}^{\mathit{K},\mathcal{H}}\right|$
Direct Form	$y\left[k\right]=g·{\sum }_0^{N-1}\left(h\left[i\right]u[k-i]\right).$	$N\mathrm{𝒶}$, $\left(N+1\right)\mathrm{𝓂}$, $\left(2N+2\right)\mathrm{𝓁}$
Direct Form Transposed	$y\left[k\right]=g·{\sum }_0^{N-1}\left(h\left[i\right]·u\right)[k-i].$	$N\mathrm{𝒶}$, $\left(N+1\right)\mathrm{𝓂}$, $\left(2N+2\right)\mathrm{𝓁}$
Symmetric	$\left\{\begin{array}{c}y\left[k\right]=g·{\sum }_0^{⌊N/2⌋}h\left[i\right]\left(u\left[i\right]+u[N-1-i]\right);\hfill \\ h\left[i\right]=h[N-1-i],i=\overline{0,⌊N/2⌋}.\hfill \end{array}\right.$	$\begin{array}{c}N\mathrm{𝒶},\left(⌊\frac{N}{2}⌋+1\right)\mathrm{𝓂},\\ \left(N+⌊\frac{N}{2}⌋+2\right)\mathrm{𝓁}\end{array}$
Anti-symmetric	$\left\{\begin{array}{c}y\left[k\right]=g·{\sum }_0^{⌊N/2⌋}h\left[i\right]\left(u\left[i\right]-u[N-1-i]\right);\hfill \\ h\left[i\right]=-h[N-1-i],i=\overline{0,⌊N/2⌋}.\hfill \end{array}\right.$	$\begin{array}{c}N\mathrm{𝒶},\left(⌊\frac{N}{2}⌋+1\right)\mathrm{𝓂},\\ \left(N+⌊\frac{N}{2}⌋+2\right)\mathrm{𝓁}\end{array}$

. The main difference between DF and DFT is that for the former, the delay word lengths are that of the input signals $u[k-i]$, while for the latter, the delays have the word length of the accumulator variable. The Symmetric and Antisymmetric cases make use of the linear phase of the filter through the regularity of the first $⌊\frac{N}{2}⌋$ coefficients as the symmetrical or antisymmetrical equivalents of the latter half of the coefficients. As mentioned in the third column of the table, this has a significant impact on the necessary multiplications and load operations involved in the implementation of $y\left[k\right],k\ge 0$. Algorithm 2 emphasizes the corresponding SISO FIR filter ISR routine starting from the mathematical basis of Equation (41) and information from Table 4.

Corollary 1.

Given a MIMO transfer matrix $\mathbf{H}\in {\mathcal{G}}_D^{p\times m}$ as in (35), where each component $H_{ij}\left(z\right)$, $i\in \overline{1,p},j\in \overline{1,m}$ can be described as an IIR filter of form (37), with second-order sections as in Table 3 or a FIR filter of form (41), with difference equations as in Table 4, the WCET can be computed as:

$${\displaystyle \mathit{WCET}}\left(\mathrm{\Psi },\mathrm{\Xi }\right)=\left(\mathcal{O}\left(1\right)+\sum _{i=1}^p\sum _{j=1}^m{\gamma }_{ij}·min\left|{\mathcal{S}\mathcal{𝓅}}^{H_{ij},\mathcal{H}}\right|\right)\times T_{clk},$$

(42)

with coefficients ${\gamma }_{ij}>0$ accounting for the hardware specifications set $\mathcal{H}$ from Definition 1, each assembly operation set ${\mathcal{S}\mathcal{𝓅}}^{H_{ij},\mathcal{H}}$ determined individually by the RCP application, given the microprocessor tick $T_{clk}>0$, and $\mathcal{O}\left(1\right)$ depends only on the other higher-priority software threads.

Algorithm 2 Finite-impulse response (FIR) filter interrupt service routine (ISR)

Input:  $H_{ij}\left(z\right)$ as in (41), topology as in Table 4 Output:  Execution time profiler for routine firFiltIsr 1: Construct software structure $H_{FIR}$ with N coefficients $h\left[i\right]$ and a Nth order delay tap 2: Initialize Nth-order input delay tap to zero 3: procedurefirFiltIsr(N,H) 4:     Read and scale $u\left[k\right]$ from input device 5:     $\tilde{u}$←$u\left[k\right]$ 6:     Shift input delay tap 7:     $\tilde{y}\leftarrow H.\mathtt{call}\left(\tilde{u}\right)$   ▷ According to column 2 of Table 4 8:     $y\left[k\right]=\tilde{y}$ 9:     Scale $y\left[k\right]$ and write to output device 10: end procedure

The proof immediately follows based on Theorem 1 by replacing the general-purpose sequence ${\mathcal{S}\mathcal{𝓅}}^{K,\mathcal{H}}$ with its corresponding sum of subsystems $H_{ij}$ from the full MIMO regulator $K\equiv \mathbf{H}$ and their definitions.

3. Case Study

The Case Study section concerns with illustrating the proposed extensions on a motor servo control example, and will encompass the following key points: process description, control performance specifications, continuous-time controller design, regulator discretization methods, sampling time optimization, selection of different controller implementation topologies along with the worst-case execution time analysis and a detailed discussion of the obtained results. The reason for selecting this example is that the proposed theory can be illustrated in the same manner for a simple process model and a modest microcontroller device setup, compared to a more complex process and an adequate microcontroller setup, leading to similar processor workloads.

3.1. Process Model and Controller Synthesis

To illustrate the proposed theoretical techniques, we consider a numerical case study based on a brushed direct-current (DC) motor position control system. Figure 5 shows the closed-loop structure, emphasizing both the process model and the 2DOF regulator structure. The motor process has a control input represented by the source voltage $V_a$ [V], along with the disturbance load torque $T_d$ [Nm], and the output is considered to be the angular position $\theta$ [rad]. As noticeable from the transfer function blocks, the DC motor model has order three, and the nominal component values are listed in

Table 5. DC motor physical parameters.

Parameter	Value	Parameter	Value
R	$2\left[\mathrm{\Omega }\right]$	L	$0.5$ [H]
$K_m$	$0.1$ [Nm·A/V2]	$K_f$	$0.2$ [Nm]
J	$0.02$ [kg · m2/s2]	$K_b$	$0.1$ [V·s/rad]

.

The process model from its two inputs to the angular position output is described by:

$$\mathrm{\Theta }\left(s\right)=\frac{1}{s}·\frac{H_L\left(s\right)}{1+H_r\left(s\right)H_L\left(s\right)H_a\left(s\right)}·T_d\left(s\right)+\frac{1}{s}·\frac{H_a\left(s\right)H_L\left(s\right)}{1+H_r\left(s\right)H_L\left(s\right)H_a\left(s\right)}·V_a\left(s\right),$$

(43)

with the armature transfer function $H_a$, along with the load component $H_L$ and the reverse loop term $H_r$ which denotes the back-electromotive voltage constant $K_b$:

$$H_a\left(s\right)=\frac{K_m}{Ls+R},H_L\left(s\right)=\frac{1}{Js+K_f},H_r\left(s\right)=K_b.$$

(44)

The 2DOF controller components have the proportional–integral–derivative plus filter (PIDF) structure for $K_{\mathrm{in}}$, while the integral term is canceled for the feedforward controller $K_{\mathrm{ff}}$. This structure allows straightforward implementation in many industrial contexts as such PIDF regulators can be directly acquired and there are multiple validated approaches in the literature for their parameter tuning [18,20]. As such, their mathematical expressions become:

$$K_{\mathrm{in}}\left(s\right)=\frac{U_{\mathrm{in}}\left(s\right)}{R\left(s\right)}=K_P+K_I\frac{1}{s}+K_D\frac{s}{T_{f}s+1},$$

(45)

$$K_{\mathrm{ff}}\left(s\right)=\frac{U_{\mathrm{ff}}\left(s\right)}{R\left(s\right)}=\tilde{b}·K_P+\tilde{c}·K_D\frac{s}{T_{f}s+1},$$

(46)

with the feedforward parameters usually specified as b and c, with $\tilde{b}=-1+b$ and $\tilde{c}=-1+c$. The command signal applicable to the motor input is comprised of two components derived from the 2DOF regulator as:

$$U\left(s\right)=U_{\mathrm{in}}\left(s\right)+U_{\mathrm{ff}}\left(s\right)=K_{\mathrm{in}}\left(s\right)E\left(s\right)+K_{\mathrm{ff}}\left(s\right)R\left(s\right).$$

(47)

The closed-loop control specifications were selected as follows: a reference tracking settling time of $t_s\le 1.5$ [s] with an overshoot $M_p\le 0.05\equiv 5$[%], with the rise time as short as possible. Additionally, regarding the disturbance rejection specifications, a load torque of 1 [Nm] must be rejected in less than $t_s^d\le 1$[s], i.e., its effect on the output measurement should become less than $0.1$ [rad] in the specified $t_s^d$, with a maximum allowed disturbance of $y_{max}^d=0.65$ [rad]. The recommended approach in such designs [12] is to tune the inner PIDF controller $K_{in}$ to account for the disturbance rejection coefficients, as $K_{ff}$ does not influence that control loop, as written in the discrete-time counterpart expressions (54). The PIDF parameter tuning has been done using global optimization methods by encompassing the desired specifications. The first iteration which covered all disturbance rejection performances was accepted and halted the optimization procedure. Additionally, a further optimization as been performed on the 2DOF parameters b and c of (46) with, again, halting the procedure after the reference tracking requirements have been fulfilled. The outcomes of the regulator tuning are illustrated in

Table 6. Regulator designs leading to the proposed 2DOF control and their corresponding coefficients.

Regulator	${\mathit{K}}_{\mathit{p}}$	${\mathit{K}}_{\mathit{i}}$	${\mathit{K}}_{\mathit{d}}$	${\mathit{T}}_{\mathit{f}}$	b	c
PID (1DOF): tracking	21.0666	8.757	7.7497	0.0014717	1	1
PID (1DOF): disturbance rejection	52.6665	70.0560	7.7497	0.0014717	1	1
PID (2DOF)	52.6665	70.0560	7.7497	0.0014717	0.4	0.2

, where, alongside the PIDF and additional PD synthesis, a separate 1DOF PIDF regulator has been also added, to account for only the tracking response.

The closed-loop step responses for the reference tracking and disturbance rejection problems are portrayed in Figure 6, showing the effects of the three controller examples from Table 6, case in which the 2DOF structure is validated, as the obtained performance metrics are $t_s\approx 1.5$[s], $M_p\approx 0[\%]$, rise time $t_r\approx 0.75$ [s] and steady-state error ${\epsilon }_{ss}=0$. As illustrated in the figure, the 2DOF structure manages to ensure both the transient response performances and the disturbance rejection behavior, being the only regulator from the proposed triplet to cover both areas.

3.2. Sampling Rate Selection

For the discretization of the PID regulators, the forward Euler method has been considered for the integrator term, with the approximation:

$$s\approx \frac{z-1}{T}=\frac{1-z^{-1}}{T{z}^{-1}},$$

(48)

while the derivative term is discretized using the backward Euler method as:

$$s\approx \frac{z-1}{Tz}=\frac{1-z^{-1}}{T},$$

(49)

leading to the expression of $K_{\mathrm{in}}$:

$$K_{\mathrm{in}}\left(z\right)=K_P+\frac{K_{I}T}{z-1}+\frac{K_D}{T_f+\frac{Tz}{z-1}},$$

(50)

with an expression of $K_{\mathrm{ff}}$ also as:

$$K_{\mathrm{ff}}\left(z\right)=\tilde{b}·K_P+\tilde{c}·\frac{K_D}{T_f+\frac{Tz}{z-1}},$$

(51)

which will be further used in the proceeding illustration of the sampling time analysis.

Starting from the continuous-time open-loop model from (43), the discrete-time equivalent using the zero-order hold method becomes:

$$\mathrm{\Theta }\left(z\right)=\mathcal{Z}\left\{{\mathcal{L}}^{-1}\left\{H_{zoh,T}\left(s\right)H_{aux}\left(s\right)\right\}\right\}·T_d\left(z\right)+\mathcal{Z}\left\{{\mathcal{L}}^{-1}\left\{H_{zoh,T}\left(s\right)H_a\left(s\right)H_{aux}\left(s\right)\right\}\right\}·V_a\left(z\right),$$

(52)

with the auxiliary notations:

$$H_{aux}\left(s\right)=\frac{1}{s}·\frac{H_L\left(s\right)}{1+H_a\left(s\right)H_L\left(s\right)H_r\left(s\right)},\mathrm{\Theta }\left(z\right)=H_{\mathrm{dist}}^{\mathrm{op}}\left(z\right)·T_d\left(z\right)+H_{\mathrm{servo}}^{\mathrm{op}}\left(z\right)·V_a\left(z\right).$$

(53)

The closed-loop system’s expression thus becomes:

$$\mathrm{\Theta }\left(z\right)=\frac{H_{\mathrm{dist}}^{\mathrm{op}}\left(z\right)}{1+H_{\mathrm{servo}}^{\mathrm{op}}\left(z\right)K_{\mathrm{in}}\left(z\right)}·T_d\left(z\right)+\frac{H_{\mathrm{servo}}^{\mathrm{op}}\left(z\right)\left(K_{\mathrm{ff}}\left(z\right)+K_{\mathrm{in}}\left(z\right)\right)}{1+H_{\mathrm{servo}}^{\mathrm{op}}\left(z\right)K_{\mathrm{in}}\left(z\right)}·R\left(z\right).$$

(54)

By imposing the functionals ${\mathcal{S}}_{K_{\mathrm{in}}}^{\overline{\mathit{\Omega }}},{\mathcal{S}}_{K_{\mathrm{ff}}}^{\overline{\mathrm{\Omega }}},{\mathcal{I}}_{K_{\mathrm{in}}},{\mathcal{I}}_{K_{\mathrm{ff}}},{\mathcal{S}}_{H_{cl}^r}^{\overline{\mathrm{\Omega }}},{\mathcal{S}}_{H_{cl}^d}^{\overline{\mathrm{\Omega }}},\mathcal{T}$, and ${\mathcal{J}}_{\mathrm{stab}\left(H_0\right)}^{{\alpha }_{\infty }}$ from Equation (23), two weighting sets $c_1$–$c_7$ were considered as in

Table 7. Functional weighting coefficients for the two considered experiments of the DC motor case study: main focus on implementability of the resulting controller versus the focus on fidelity compared to the continuous-time control counterpart.

Parameter	Implementability Case Value	Fidelity Case Value
$c_1$	0.1	100
$c_2$	0.1	100
$c_3$	2000	1
$c_4$	2000	1
$c_5$	0.1	300
$c_6$	0.1	100
$c_7$	200	1
$T_s^{opt}$	$T_{s,1}^{opt}=2.866\times {10}^{-3}$ [s]	$T_{s,2}^{opt}=1.260\times {10}^{-4}$ [s]

, corresponding to two distinct experiments, the first numerical column designating emphasis on the difficulty of implementation functionals, while the second numerical column coefficients focus mainly on open-loop and closed-loop fidelity. With said coefficients, using a general-purpose PSO implementation as specified in Remark 3 [28], the optimal implementability sampling period becomes $T_{s,1}^{opt}=2.866\times {10}^{-3}$ [s], while the fidelity sampling period is obtained at $T_{s,2}^{opt}=1.260\times {10}^{-4}\left[\mathrm{s}\right]\ll T_{s,1}^{opt}$. To further extend the analysis, two extra sampling periods will be added in the discussion, the first representing the value obtained by applying the classical Shannon–Nyquist theorem (9), leading to a sampling period smaller than half the least time constant of the regulator, i.e., $T_{s,3}<\frac{T_f}{2}\Rightarrow T_{s,3}=7.0081\times {10}^{-4}$ [s], while the latter relevant value is considered to be $T_{s,4}=T_{s,1}·\left(1+0.03\right)=2.9520\times {10}^{-3}$, which represents a $3[\%]$ disturbance increase on the ideal implementability sampling rate, which causes instability in the closed-loop system. Figure 7 gathers all functionals and their weighted sums in its six subfigures, while also marking the positions of $T_{s,1}^{opt}$, $T_{s,2}^{opt}$ and $T_{s,3}$. The open and closed-loop similarity functionals $J_1$, $J_2$, $J_5$, $J_6$ are non-monotonic with respect to a variable $T_s$, while the implementability functionals $J_3$, $J_4$, $J_7$ are principally monotonically decreasing. The existing exceptions appear due to numerical errors.

There are multiple approaches for implementing the two PID regulators. Both $K_{\mathrm{in}}\left(z\right)$ and $K_{\mathrm{ff}}\left(z\right)$ can be fully encompassed into a biquadratic filter topology and, furthermore, a first-order structure for $K_{\mathrm{ff}}$. The two main ones, given the simplicity of their structure, would be to implement it in parallel versus in series. For the inner regulator, the parallel topology, denoted with the superscript p can be split into three subsystems:

$$K_{\mathrm{in}}^p\left(z\right)=H_{in,1}\left(z\right)+H_{in,2}\left(z\right)+H_{in,3}\left(z\right)=K_P+\frac{K_{I}T}{z-1}+\frac{K_D(z-1)}{(T_f+T)z-T_f},$$

(55)

while the series topology, denoted with the superscript s, is:

$$K_{\mathrm{in}}^s\left(z\right)=\frac{b_2{z}^2+b_{1}z+b_0}{z^2+a_{1}z+a_0}=g·\frac{\widetilde{b_2}{z}^2+{\tilde{b}}_{1}z+{\tilde{b}}_0}{a_2{z}^2+a_{1}z+a_0},a_2=1.$$

(56)

where the equivalent normalized coefficients are obtained:

$$b_2=\frac{K_P(T_f+T)+K_D}{T_f+T};b_1=\frac{-K_P(2T_f+T)+K_{I}T(T_f+T)-2K_D}{T_f+T};$$

(57)

$$b_0=\frac{K_P{T}_f-K_{I}T{T}_f+K_D}{T_f+T};a_1=-\frac{2T_f+T}{T_f+T};a_0=\frac{T_f}{T_f+T}.$$

(58)

In the same manner, the parallel form of the feedforward regulator is adapted as:

$$K_{\mathrm{ff}}^p\left(z\right)=\tilde{b}{H}_{in,1}\left(z\right)+\tilde{c}{H}_{in,3}\left(z\right)=\tilde{b}·K_P+\tilde{c}·\frac{K_D(z-1)}{(T_f+T)z-T_f},$$

(59)

with a first-order series form of:

$$K_{\mathrm{ff}}^s\left(z\right)=\frac{b_{1}z+b_0}{z+a_0}=g·\frac{{\tilde{b}}_{1}z+{\tilde{b}}_0}{a_{1}z+a_0},a_1=1.$$

(60)

and equivalent normalized coefficients:

$$b_1=\frac{\tilde{b}{K}_P(T_f+T)+\tilde{c}{K}_D}{T_f+T};b_0=\frac{-\tilde{b}{K}_P{T}_f-\tilde{c}{K}_D}{T_f+T};a_0=\frac{-T_f}{T_f+T}.$$

(61)

3.3. Execution Time Analysis

After the discretization procedure as specified by (48) and (49), using the right-hand side notation for the coefficients in (56) and (60), the numerical values of the coefficients become as in

Table 8. DC motor case study discrete-time ideal coefficients using the deduced sampling rates.

Coefficient	${\mathit{T}}_{\mathit{s},1}^{\mathbf{opt}}$	${\mathit{T}}_{\mathit{s},2}^{\mathbf{opt}}$	${\mathit{T}}_{\mathit{s},3}$	${\mathit{T}}_{\mathit{s},4}$
g: $K_{\mathrm{in}}\left(z\right)$	$5.3184815\times {10}^3$	$5.3184815\times {10}^3$	$5.3184815\times {10}^3$	$5.3184815\times {10}^3$
${\tilde{b}}_2$: $K_{\mathrm{in}}\left(z\right)$	1	1	1	1
${\tilde{b}}_1$: $K_{\mathrm{in}}\left(z\right)$	$-1.980677$	$-1.999150$	$-1.995275$	$-1.980097$
${\tilde{b}}_0$: $K_{\mathrm{in}}\left(z\right)$	$0.980751$	$0.999150$	$0.995279$	$0.980175$
$a_1$: $K_{\mathrm{in}}\left(z\right)$	$-0.052546$	$-1.914358$	$-1.523809$	$0.005877$
$a_0$: $K_{\mathrm{in}}\left(z\right)$	$-0.947453$	$0.914358$	$0.523809$	$-1.005877$
g: $K_{\mathrm{ff}}\left(z\right)$	$-4.244251\times {10}^3$	$-4.244251\times {10}^3$	$-4.244251\times {10}^3$	$-4.244251\times {10}^3$
${\tilde{b}}_1$: $K_{\mathrm{ff}}\left(z\right)$	1	1	1	1
${\tilde{b}}_0$: $K_{\mathrm{ff}}\left(z\right)$	$-0.985500$	$-0.999362$	$-0.996454$	$-0.985065$
$a_0$: $K_{\mathrm{ff}}\left(z\right)$	$0.947453$	$-0.914358$	$-0.523809$	$1.005877$

, using the for deduced sampling rates from Section 3.2, $T_{s,1}^{opt}$, $T_{s,2}^{opt}$, $T_{s,3}$, $T_{s,4}$. As observed, the gain coefficient $g\equiv b_2$ for $K_{\mathrm{in}}$ and $g\equiv b_1$ for $K_{\mathrm{ff}}$, respectively, remains invariant in this set of experiments, while the remaining non-unitary coefficients vary with respect to $T_s>0$ and necessitate increasingly more decimals as the sampling rate tends to zero.

To account for the necessary word length analysis, two frequently used configurations have been considered: the first case is to store the operands, i.e., coefficients and inputs, states, outputs into 16-bit registers, considered the standard length for the RISC machine hosting the 2DOF controller, followed by a set of 32-bit length registers, which will increase the working precision, but with added execution time overhead, as modeled in continuation. Given the dynamic range of the final filter gains g, this final multiplication will be considered separately. Thus, the set $\mathcal{H}$ for this case study will encompass the properties:

• it has a clock tick period of $T_{clk}=\frac{1}{1\left[\mathrm{MHz}\right]}=1$ [$\mathsf{\mu }$s], obtainable on standard microcontrollers by configurating the phase-locked loop circuit to a lower-power setting;
• implements the second-order sections using the Direct Form II for $K_{\mathrm{in}}$ and a series connection for the first-order term $K_{\mathrm{ff}}$, denoted $y\left[k\right]=f\left(y[k-i],u[k-i]\right)$;
• two configurations for the word lengths of the operands: 16-bit and 32-bit;
• applies saturation on output command signals $y\left[k\right]$;
• applies anti-windup on the integrator term of $K_{\mathrm{in}}$ using the back-calculation method:

  $${\mathcal{w}}_{\underline{I},\overline{I}}\left(x[k+1]\right)=x\left[k\right]+\left[K_w·\left({\mathcal{s}}_{\underline{Y},\overline{Y}}\left(y\left[k\right]\right)-y\left[k\right]\right)\right]·\left(K_{i}T\right)·u\left[k\right],$$

  (62)

  with the additional parameters $K_w$, $\underline{I}$, $\overline{I}$ represented using a 16-bit word length;
• the output measurement $\mathrm{\Theta }$ is gathered as a sum of impulses, with a maximum expected frequency of $f_{max}=$ 80,000 impulses per 100 [ms] time unit, and each such impulse triggers a hardware interrupt with a 15 assembly operation stack commutation cost; the scaling is then performed using a multiplication with a 16-bit variable;
• $K_{\mathrm{ff}}$ accepts the reference signal $r\left[k\right]$ as input, while $K_{\mathrm{in}}$ accepts the error $e\left[k\right]=r\left[k\right]-\mathrm{\Theta }\left[k\right]$.

The mathematical operations $p_i\in \mathcal{O}$, along with their assembly instruction correspondents as in (29) are detailed in

Table 9. Encountered operations in the implementation of the 2DOF structure for the DC motor case study in the hypothesis of a base word length of $L=16$ bits and two considered quantization levels for the controller coefficients: 16-bit and 32-bit lengths, respectively.

${\mathit{p}}_{\mathit{i}}=\left(\#{\mathit{r}}_{\mathit{i}}\right)$	${\mathit{p}}_{\mathit{i}}$	${\mathit{K}}_{\mathbf{in}}^{\mathit{s}}$	${\mathit{K}}_{\mathbf{in}}^{\mathit{p}}$	${\mathit{K}}_{\mathbf{ff}}^{\mathit{s}}$	${\mathit{K}}_{\mathbf{ff}}^{\mathit{p}}$	$\begin{array}{c}{\mathit{\gamma }}_{\mathit{i}}\\ \mathit{16}\mathit{b}\end{array}$	$\begin{array}{c}{\mathit{\gamma }}_{\mathit{i}}\\ \mathit{32}\mathit{b}\end{array}$	Observations for$\mathit{K}\left(\mathit{z}\right),\mathcal{H}$
1	$\mathrm{𝓁}$	14	18	8	13	1	2	$y\left[k\right]=f\left(y[k-i],u[k-i]\right)$
1	$\mathrm{𝒶}$	4	5	2	3	1	2	$y\left[k\right]=f\left(y[k-i],u[k-i]\right)$
1	$\mathrm{𝓂}$	6	5	4	4	1	6	$y\left[k\right]=f\left(y[k-i],u[k-i]\right)$
8	$\mathrm{𝓈}$	1	1	1	1	1	2	$\mathrm{𝓈}\left(y\left[k\right]\right)$
20	$\mathrm{𝓌}$	1	1	0	0	1	6	Integrator anti-windup
$120\times {10}^3\times T_s$	$\mathrm{𝒶}$	1	1	0	0	1	1	$\mathrm{\Theta }\left[k\right]$ impulse counter
1	$\mathrm{𝓂}$	1	1	0	0	1	6	$\mathrm{\Theta }\left[k\right]$ scaling
1	$\mathrm{𝒶}$	1	1	0	0	1	2	$e\left[k\right]=r\left[k\right]-\mathrm{\Theta }\left[k\right]$
1	$\mathrm{𝓁}$	3	3	0	0	1	2	$e\left[k\right]=r\left[k\right]-\mathrm{\Theta }\left[k\right]$
1	$\mathrm{𝓁}$	2	2	1	1	1	2	$y[k-i]\leftarrow y[k-i-1]$
1	$\mathrm{𝓁}$	2	2	1	1	1	2	$u[k-i]\leftarrow u[k-i-1]$

, with an emphasis on each type of operation based on its physical significance, as written in the last column. The hyperparameters ${\gamma }_i$ in the case of standard 16-bit word lengths have been considered with the value 1, denoting that each base arithmetical operations costs only $T_{clk}$, while the values are scaled upwards for the case when the operands exceed the word length to 32 bits. The impulse counter assembly operation cost for $\mathrm{\Theta }\left[k\right]$ has been computed as 80,000 $\times 0.1\times 15\times T$. Additionally,

Table 10. DC motor case study WCET analysis in the conditions from Table 9, by emphasizing the number of assembly operations $\left|{\mathcal{S}\mathit{𝓅}}^{K,\mathcal{H}}\right|$ necessary to implement the 2DOF structure comprised of $K_{\mathrm{in}}$ and $K_{\mathrm{ff}}$, along with the processor usage level for the three stable sampling rate values.

	$\begin{array}{c}\mathit{16}\mathit{b}\\ {\mathit{K}}_{\mathbf{in}}^{\mathit{s}}\end{array}$	$\begin{array}{c}\mathit{16}\mathit{b}\\ {\mathit{K}}_{\mathbf{ff}}^{\mathit{s}}\end{array}$	$\begin{array}{c}\mathit{16}\mathit{b}\\ {\mathit{K}}_{\mathbf{in}}^{\mathit{p}}\end{array}$	$\begin{array}{c}\mathit{16}\mathit{b}\\ {\mathit{K}}_{\mathbf{ff}}^{\mathit{p}}\end{array}$	$\begin{array}{c}\mathit{32}\mathit{b}\\ {\mathit{K}}_{\mathbf{in}}^{\mathit{s}}\end{array}$	$\begin{array}{c}\mathit{32}\mathit{b}\\ {\mathit{K}}_{\mathbf{ff}}^{\mathit{s}}\end{array}$	$\begin{array}{c}\mathit{32}\mathit{b}\\ {\mathit{K}}_{\mathbf{in}}^{\mathit{p}}\end{array}$	$\begin{array}{c}\mathit{32}\mathit{b}\\ {\mathit{K}}_{\mathbf{ff}}^{\mathit{p}}\end{array}$
$\left|{\mathcal{S}\mathit{𝓅}}^{K,\mathcal{H}}\right|$	77	24	81	30	246	64	250	76
$\left|{\mathcal{S}\mathit{𝓅}}^{\left\{K_{\mathrm{in}},K_{\mathrm{ff}}\right\},\mathcal{H}}\right|$	101		111		310		326
WCET$T_{s,1}^{opt}=2866$[$\mathsf{\mu }$s]	429[$\mathsf{\mu }$s]		439[$\mathsf{\mu }$s]		638[$\mathsf{\mu }$s]		654[$\mathsf{\mu }$s]
WCET$T_{s,2}^{opt}=126$[$\mathsf{\mu }$s]	101[$\mathsf{\mu }$s]		111[$\mathsf{\mu }$s]		310[$\mathsf{\mu }$s]		326[$\mathsf{\mu }$s]
WCET$T_{s,3}=700.8$[$\mathsf{\mu }$s]	170[$\mathsf{\mu }$s]		180[$\mathsf{\mu }$s]		379[$\mathsf{\mu }$s]		395[$\mathsf{\mu }$s]
$\mathcal{U}\left(\mathrm{\Psi },\mathrm{\Xi },T_{s,1}^{opt}\right)$	$14.96[\%]$		$15.31[\%]$		$22.26[\%]$		$22.81[\%]$
$\mathcal{U}\left(\mathrm{\Psi },\mathrm{\Xi },T_{s,2}^{opt}\right)$	$80.15[\%]$		$88.09[\%]$		$246.03[\%]$		$258.73[\%]$
$\mathcal{U}\left(\mathrm{\Psi },\mathrm{\Xi },T_{s,3}^{opt}\right)$	$24.25[\%]$		$25.68[\%]$		$54.08[\%]$		$56.36[\%]$

totalizes the number of assembly operations based on the previous table’s description, along with computing the worst-case execution times for the three stable sampling rate values: $T_{s,1}^{opt}$–$T_{s,3}$. As seen from the table, the optimal sampling rate solution featuring implementability emphasis $T_{s,1}^{opt}$ occupies the microprocessor for less than $25[\%]$ of its capability in either 16 or 32-bit quantizations alike, with the Shannon theorem approach $T_{s,3}$ following with sufficient headroom in the scheduler algorithm, while the fidelity-based approach in this case occupies the scheduler with small margins in the case of 16-bit quantizations, and exceeds the allowed time frame in the case of the 32-bit configuration.

Based exclusively on the previous WCET analysis, there are several feasible solutions. In order to decide between the three sampling rate and quantization pair configurations, further analysis is performed on the frequency response of the $K_{\mathrm{in}}$ and $K_{\mathrm{ff}}$ controllers, along with the closed-loop responses using said controllers. As such, Figure 8 and Figure 9, respectively, gather the previously said behaviors. In addition to the three sampling rates, for completeness, the fourth, unstable sampling rate $T_{s,4}$ is shared, along with the ideal continuous-time equivalent controllers in the frequency-response plot. In both figures, columns one and two expose the 32-bit and 16-bit quantization configurations, respectively, while the lines distinguish between the $K_{\mathrm{in}}$, $K_{\mathrm{ff}}$ controllers in the frequency-response figure and reference step response compared to the step disturbance responses in the time-domain simulation, respectively. The main conclusion drawn from this final pair of Figures is that the considered 2DOF control scheme is sensitive to the coefficient quantization levels, such that for $T_{s,2}$ and $T_{s,3}$, the only acceptable solution would be to use the precise 32-bit configuration, with the exception of the implementability solution which manages to follow the imposed closed-loop transient response specifications with a reasonable degradation of the performances. The 16-bit quantization frequency responses for $K_{\mathrm{in}}\left(z\right)$ drastically alter the integrator effect, which, in effect, disturbs the ability of the closed-loop motor system to track reference signals and to reject step-like load torque disturbances.

To conclude the case study, the acceptable solutions are to consider the implementability-based sampling rate optimum $T_{s,1}^{opt}$, with practically ideal behavior if a 32-bit word length setup is acceptable, and with a small performance degradation if only the 16-bit standard word length is supported, depending also on other execution threads running in the microprocessor, not covered in this experiment.

Further extensions, as in using more complex controllers for both the tracking regulator and the feedforward component, considering cases such as fractional-order controllers or 2DOF robust control synthesis results can be treated in an analogous manner by splitting such control laws into their component second-order sections and applying the theory from Table 3 and Table 4, along with Algorithms 1 and 2, respectively.

4. Conclusions

This paper gathered a set of analysis and design tools to determine the sampling rate of one and 2DOF control structures using an optimization-based approach, along with an approach of deducing a WCET for linear and time-invariant-based regulators through a formal language model which can be implemented in an RCP software tool. The execution time model is based on a deterministic, RISC architecture, where each operation is quantified with the same base clock tick duration. The end-to-end DC motor case study emphasizes the design of the controllers for the widely used benchmark system, by also focusing on the proposed framework.

Future work will concern on an online sampling rate optimization, as it is necessary in linear-parameter varying and linear-time variant models, along with the study and analysis of the continuity property of the regulator coefficient quantization effects as a function of the sampling rate. The problem of process controllability loss is also known in the literature and will be investigated for variable sampling rates in subsequent research. The mathematical framework proposed and extended in this paper will be included in the software toolbox initially proposed by the authors in [22], with the great advantage of obtaining an end-to-end solution for RCP, starting from the continuous-time controller design up to the optimal numerical implementation of said controller on a microprocessor-based system with a given set of specifications. A separate theoretical direction would be to investigate the link between sampling rate $T>0$ and quantization step $Q>0$ to guarantee that the minimum imposed performances are also fulfilled in the discrete domain.