Interpolation-based, minimum-time piecewise constant control of linear continuous-time SISO systems

We present a novel approach to the system inversion problem for linear, scalar (i.e. single-input, single-output, or SISO) plants. The problem is formulated as a constrained optimization program, whose objective function is the transition time between the initial and the final values of the system’s output, and the constraints are (i) a threshold on the input intensity and (ii) the requirement that the system’s output interpolates a given set of points. The system’s input is assumed to be a piecewise constant signal. It is formally proved that, in this frame, the input intensity is a decreasing function of the transition time. This result lets us to propose an algorithm that, by a bisection search, finds the optimal transition time for the given constraints. The algorithm is purely algebraic, and it does not require the system to be minimum phase or nonhyperbolic. It can deal with time-varying systems too, although in this case it has to be viewed as a heuristic technique, and it can be used as well in a model-free approach. Numerical simulations are reported that illustrate its performance. Finally, an application to a mobile robotics problem is presented, where, using a linearizing pre-controller, we show that the proposed approach can be applied also to nonlinear problems.


Introduction
Systems inversion, or output tracking, is a long-studied topic in control theory since the seminal work of Brockett and Mesarović [1], where the term reproducibility was used to denote it. Broadly speaking, it consists of finding, for a given system, the input that produces the desired output (thus justifying the name system inversion).
This technique is most often applied to systems already equipped with a feedback controller, which cares for stability and disturbance rejection, but need a better tracking in the transient phase, that makes a feed-forward action appropriate. The systems inversion problem has been solved for linear SISO systems by Brockett [2], and since then it has been generalized to multivariable linear systems by several authors, e.g. Sain and Massey [3], and Silverman and Payne [4].
Motivated by practical problems in the control of flexible structures, the extension to multivariable nonlinear systems has been realized: see Hirschom's paper [5] for the theoretical aspects, and Di Benedetto and Lucibello's [6] for an application to the endpoint control of a flexible-arm manipulator. In the latter, the system under consideration is also time-varying.
More recently, the system inversion problem has been cast into an optimization frame; for example, Iamratanahl and Devasia [7] have investigated the minimum-time and minimum-energy output tracking problem for linear, multiinput, multi-output (MIMO) systems, yet constrained to square plants, i.e. plants having the same number of inputs and outputs. Afterwards, in [8] it has been shown that the output minimum transition time for nonlinear square MIMO systems under a constraint on the input intensity can be improved beyond the performance of the 'bang-bang control' by using pre-and post-actuation.
On returning to linear SISO systems, which are the subject of this paper, the system inversion problem has been often cast into transferring the system's output from an initial value to a final one with some additional criteria, such as constraints on the shape of the output (for instance, not peeking, or having bounded derivatives [9]), or using minimum input energy [10].
The approach of 'imposing the shape' of the output taken in [9] has been extended recently in [11], where nonlinear scalar systems are considered, with the constraints fulfilment embedded into the input design optimization problem, and in [12], where the behavioural approach has been exploited to minimize the transition time and enforce a given smoothness degree of the input.
The above problems are made more difficult by the presence of right-half-plane zeros in the plant, which are common, for example, in the linearized models of flexible structures and cannot be recovered by the feedback controller: see [13] for an application. As explained in that paper, the transient tracking performance of linear MIMO systems is limited by its right-half-plane zeros, and a measure of performance that involves those zeros only can be given.
The limit case of nonhyperbolic systems, i.e. systems having zeros on the imaginary axis, which can lead the input signal to be very large [14], has been considered by Jetto et al.: see [15,16]. In particular, the latter is one of the very few papers we are aware of that casts the inversion problem into a discrete-time frame, and assumes an 'a priori' structure of the input as a spline function. See also Marconi et al. [17] for a geometric approach to the problem for linear, discrete-time SISO systems.
In this paper, for SISO time-invariant systems we present an interpolation approach to the output tracking problem: the system's output is forced to pass through a set of points taken from an ideal output, using a piecewise-constant input that minimizes the transition time, subject to a constraint on its intensity.
As already said, constraining the input intensity is often considered in the inversion literature, owing to its relevance in practical problems: see, e.g. [18], where a saturation approach is used, or [19], where constraints on both the input and its change rate are imposed.
The use of piecewise-constant input functions, also known as block-pulse functions (BPF), gives some advantages [20], or obeys to some constraints [21], where the structure of the controller is concerned. Piecewise-constant input functions are widely used, for example, in mobile robotics, where robots usually travel for the most part of their journey at a constant speed.
The main result of this contribution is a locally exact relationship between the transition time and the control effort needed for the transition. That relationship allows us to propose an algorithm that provides the minimum transition time for a given number of output values to be interpolated and a given bound on the input, and computes the values of the control signal segments.
Unlike the previously cited papers, owing to its purely algebraic nature, the algorithm we propose does not rely upon the stability and/or the minimum-phase property of the controlled system. It can be applied to time-varying systems as well, although in that case the relationship between the transition time and the control effort does not apply rigorously, and it has to be viewed as a heuristic approach. It is also worth to remark that the algorithm can be used in a modelfree context, because it only requires the knowledge of the impulse response of the sampled data system: see Remark 1 in the next section.
A discussion about the behaviour of the system's output far from the interpolation points is given, and a numerical comparison of the algorithm performance versus a method from the recent literature is also presented, which shows the effectiveness of the proposed approach.
Finally, an application of the proposed control approach to mobile robotics is presented, where a unicycle robot is driven to track a given trajectory. The method is applied to the input-output linearized model obtained by a feedback linearization technique [22], and then the optimal signal is used to feed the nonlinear system. This suggests that our approach can be straightforwardly applied also to nonlinear inversion problems when a linearizing-decoupling pre-controller is available.

Problem statement
Let be the continuous-time, linear time-invariant SISO system of order n described by the state-space model where x ∈ R n is the state vector, u ∈ R and y ∈ R are the input and output signals, respectively, with A ∈ R n×n , b and c ∈ R n . We assume x(0) = 0 as the initial state, which in turn implies y(0) = 0. For a given T > 0, let t k = kT , k ∈ N, be a set of equally spaced time instants, and assume that the input u is constant in each interval [t k , t k+1 ), i.e. u(t) ≡ u(t k ) in that interval. Then the output of system at the time instants t k > 0 is: Given a set {y * k , k = 1, . . . , m} of desired output values, find the smallest interval size T and the corresponding input sequence {u(t 0 ), . . . , u(t m−1 )}, such that y(t k ) = y * k , and a threshold on the input intensity is fulfilled: In other words, the interval size T * has to be chosen that solves the optimization problem: where y(t k ), k = 1, . . . , m, are given by Eq. (1).
To solve the problem we first note that, as the matrix A is constant, we have, independently of k: and, denoting by θ T the integral at the r.h.s. of the previous formula, we obtain by iteration: The sequence {h p (T ), p ∈ N}, which depends on the choice of T , is the impulse response sequence of the sampled-data system corresponding to the continuous time system for that value of T . Also notice that, if we let

Remark 1
The generic value h p (T ) can be seen as the difference between two consecutive values of the step response y s of the system, i.e.: This allows to obtain the values h p (T ) also by a real-ground experiment or a numerical simulation.
In view of Eq. (3), the constraint y(t k ) = y * k is: which is a linear algebraic system of the form where H T is the lower-triangular Toeplitz's matrix If h 1 (T ) = 0, then the system H T U = Y * has a unique solution, and we rewrite Problem (2) in the more compact form

Main result
It is intuitive, on the basis of physical considerations, that the smaller we take T , the larger the maximum input size U ∞ becomes; this norm is bounded by [23, p. 345] Where the infinity norm of H −1 T is concerned, we will prove the following result.

Proposition 1 For sufficiently small values of T , the infinity norm of H −1 T obeys to
where q is the order of first derivative of h that does not vanish at the origin, and η is independent of T .
Proof See Appendix.

Remark 3
It is worth noticing that the value q in Proposition 1 is equal to the relative degree of system.
As a consequence of this result and of Eq. (8), increasing T makes, as expected, the upper bound on the control effort smaller, thus enlarging the feasible region of Problem (7), provided that T is small enough for Proposition 1 to hold. This is the rationale behind the algorithm to solve Problem (7) presented in the next section.
The numerical solution of Toeplitz's linear system H T U = Y * can be carried out either by any standard technique for triangular systems, with a computational cost of O(m 2 ) elementary operations, or by the tailored Commenges-Monsion's algorithm [24], which requires O(m log(m)) flops.

Algorithm
The result of Sect. 3 enables us to propose an iterative algorithm to solve Problem (7). The idea is to start the search for the solution from a small value of T , that will probably violate the constraint on the input, and then increase that value by repeatedly doubling it until the constraint is satisfied.
In view of Proposition 1, each time T is doubled the intensity of the input needed for the control decreases, until a feasible solution is found, provided that T is small enough for Proposition 1 to apply. Formally, we have the following result about the local convergence of the algorithm. (7) is small enough for Proposition 1 to hold, owing to the continuity of the objective function, the sequence generated by the proposed algorithm converges to it.

Algorithm 1 Bisection-type search
The above result, however, cannot guarantee the convergence to the solution of any instance of Problem (7), because the optimal solution could exist, but be large enough that the approximation given by Proposition 1 does not hold (that essentially depends on the time constants of the system). In that case the proposed algorithm has to be viewed as a heuristic procedure.

Output behaviour analysis
A limitation of the proposed method is that it guarantees the exact matching between the real output and the ideal one at the time instants t k , but not in the inner points of any interval [t k , t k+1 ), because no specification is given for those points.
However, it is generally desirable that the values of the output in the open interval (t k , t k+1 ) are not too far from the values y * k attained at t = t k , or similarly at t = t k+1 . It is intuitive, and it has been shown in [25] that this can be obtained by making T small enough, provided that this is allowed by the constraints. Here we give a limitation on the distance which accounts for the variability of the output to the right of each interpolation point.
In that interval, the system's output is given by and taking its first-order Taylor approximation around the instants t k we get, for t ∈ [t k , t k+1 ), The computation ofẏ(t k ) yieldṡ and, by taking into account the constraint on the input size, we get which is a bound that can be computed independently of the goal Y * to be reached (but it needs that the interval length T has been fixed), or which is a sharper bound that can be obtained once the target Y * has been fixed and the input has been computed. Denoting by Q 1 (k, T , U M ) the first sum and by Q 2 (k, T , U ) the second one, we have two bounds for e k (t): which are clearly null in t = t k , and maximized by t → t k+1 ; in the latter case we have Thus, the first formula allows, for a given T and a given threshold U M on the input size, to know the worst a priori error within each interval.

The time-varying case
The proposed approach can be extended to time-varying linear SISO systems described by with the initial condition x(t 0 ) = 0. Given again a grid of equally spaced points t k = kT , k ∈ N, and a piecewise constant input sequence {u(t 0 ), . . . , u(t m−1 )}, we can write where is the state transition matrix [26] of the system. Denoting by M(t k , t k−1 ) the integral in the above formula we have and we can again iterate, starting from x(t 0 ) = 0, to obtain the solution at the time instants t k where the composition property (t, σ ) (σ, ξ ) = (t, ξ) and (t, t) = I have been used [26].
Defining t i ), the interpolation condition can be written also for the system t which is again a lower triangular system, although the matrix is not a Toeplitz's anymore. This implies that the proposed algorithm can still be used to find the optimal input sequence, but it has to be regarded as a heuristic, because the result of Proposition 1 does not apply.

Numerical experiments
In this section we report the results of three numerical experiments, aimed at investigating the effectiveness of the proposed algorithm in solving the minimum-time interpolation problem. The computations have been performed in MATLAB running on a standard PC. We denote by τ = mT the length of the control interval, and in all experiments we take T min = 10 −8 and ε = 10 −6 .
In the three experiments, the ideal output is the polynomial function defined in the interval [0, τ ] proposed in [9] y d (t, τ ) = That function and its first p derivatives vanish in t = 0, the first p derivatives vanish also in t = τ , and the function is monotonically increasing independently of p, thus representing a smooth transition from the initial to the final values of the output. The parameters to be set for each experiment are the input threshold U M , the number of desired output values m and the parameter p of (12). Following [9], the latter will be taken equal to the relative degree of the plant.

Remark 4
It is worth reporting that the three experiments have been performed also using the model-free approach illustrated in Remark 1, obtaining exactly the same results shown in this section.

Experiment 1
In the first numerical experiment we consider the problem proposed in [27], where a third-order nonminimum phase the ideal output function y d (t, τ ) from which we have taken the desired output values y * k is polynomial (12) with p = 1, i.e.
The approach pursuit in that paper relies on the concept of cause/effect pairs associated with the system, and it is aimed at finding, between all input/output pairs satisfying some smoothness property, the one that minimizes the transition time between two set points. In this experiment we set m = 8, U M = 1.4. Figure 1 shows the error y(t) − y d (t, τ ) and the piecewise constant control signal. The minimum control time found by the proposed algorithm is τ * = 0.47 s, which is half the value τ * = 0.94 s obtained in [27] using U M = 1.2.
We remark that a smooth control signal is generated by the method in [27], differently from our approach that returns a BPF, which is less expensive to be generated.

Experiment 2
In the second experiment we have considered the problem proposed in [9], where a plant with transfer function P(s) = 377 (s + 2) ((s + 2) 2 + 9)((s + 3) 2 + 49) with K c = 3.6172, δ = 0.4323 and ω n = 5.1073, in a unityfeedback scheme that gives rise to a sixth-order system with a relative degree of three. The value p = 3, equal to the relative degree of the closedloop plant, has been taken in Eq. (12), yielding to the ideal output function This problem has been solved in [9] using a bisection search scheme similar to the one proposed here, which involves a global optimization search in a single variable at each step, and there the optimal value τ * = 0.367 s has been found, using a threshold U M = 3 on the input intensity.
We have used for our algorithm the threshold value U M = 2.8 and m = 5, and have obtained τ * = 0.358 s, which is slightly better than the value obtained in the cited paper. Figure 2 shows the functions y(t) and y d (t, τ ) and the piecewise constant control signal. There you can see that, owing to a couple of complex poles with very low damping (δ = 0.15), the system output peeks in the time intervals between the interpolation points.
The last comment made for the previous numerical experiment also applies to this one.

Experiment 3
In the third experiment we consider the problem presented in [28], where a rocket must follow a one-dimensional position profile, and the rocket's mass is time-varying because of the fuel consumption during the journey.
It is assumed that the rocket's mass varies according to where m R is the mass of the empty rocket, m F is the initial mass of the fuel, and α is the fuel's consumption rate. On denoting by z(t) the position of the rocket and by ζ the viscous friction coefficient, the rocket's motion is governed and, returning to the dot notation for the time derivatives: Suppose that a proportional position controller with value K has been pre-wrapped to the plant, obtaining the model
The plant's physical parameters are m R = 1.5 kg, m F = 1 kg, ζ = 10 −5 N s/m, α = 0.5 s −1 , and the proportional controller is K = 5. The ideal output trajectory in this case is given again by Eq. (12) with p = 4, which is an arbitrary choice, although the concept or relative degree has been introduced also for time-varying systems [30]. The numerical computation of functions h i j of Eq. (11) for this model has been performed by a simulation, using the integrator ode45 of MATLAB, assuming as input a suitable set of m block-pulse functions. The parameters of the algorithm we have used in this case are m = 8 and U M = 12. The optimal value of τ found by our method is τ * = 1.41 s, and the tracking error and the input are shown in Fig. 3.

An application to mobile robotics
In this section, the application of the proposed approach to the control problem of the unicycle mobile robot is considered. That kind of mobile robot is modelled as in [22] ⎧ ⎨ ⎩ẋ where x r , y r , θ r represent the robot's pose, i.e. the robot centre's coordinates in the plane and its heading, while v r and ω r are the robot's inputs. In particular, v r is the driving velocity, i.e. the modulus (including the sign) of the velocity vector of the contact point, whereas the steering velocity ω r is the angular velocity of the wheel around the vertical Fig. 5 The proposed control scheme for the unicycle axis passing through the robot's centre: see Fig. 4. Assuming that the robot's inputs are constrained by |v r (t)| ≤ V and |ω r (t)| ≤ , the goal is to design a control law that drives the robot along a desired trajectory (x * (t), y * (t)).
To this aim, the standard input-output linearization method [22] has been used: by applying the control law robot's model (14) is straightforwardly reformulated as where (x B , y B ) are the coordinates of a point B located along the sagittal axis of the unicycle at a distance b > 0 from its centre. Note that the robot's heading angle θ r cannot be controlled by this technique. The proposed method is then applied to linearized and decoupled model (16) to drive the robot along a desired path; the resulting control scheme is shown in Fig. 5. The constraints to impose on the 'virtual inputs' u x and u y to satisfy the real constraints on v r and ω r can be obtained by writing linearization law (15) explicitly v r (t) = u x (t) cos(θ r (t)) + u y (t) sin(θ r (t)) ω r (t) = (u y (t) cos(θ r (t)) − u x (t) sin(θ r (t))) b −1 and obtaining It follows easily that by imposing the 'virtual constraints' |u x (t)| ≤ U M and |u y (t)| ≤ U M , with U M ≤ min{V , b}/2, also the real constraints will be satisfied.
The BPF control laws are then designed independently on the x B and y B axes using single integrator-based evolution (16). Two candidate values for T are obtained, namely T x and T y , one for each axis. To be compliant with the desired trajectory on both axes, the value of T is chosen as the largest of them, and then the virtual signals u x (t) and u y (t) are computed. The resulting control inputs u x (t) and u y (t) are reported in Fig. 7, and as expected, they are compliant with the imposed constraints.
The evolution of (x B , y B ) coincides with the reference trajectory at the points (x * (kT ), y * (kT )), while the robot's centre trajectory error satisfies as it is expected due to the input-output linearization.

Conclusions
We have proposed a control methodology to solve the inversion problem for linear SISO systems with a minimum-time criterion and a bounded piecewise-constant input. We have shown that, in the proposed setting, the norm of the input is a decreasing function of the control interval length, a property that is often assumed true on the basis of physical considerations but that is formally proven only in specific cases.
Two distinguishing features of the proposed method are that it does not rely on any assumption, besides linearity, on the plant, and that it can be used in a model-free context, as the only data needed are the step response of the system.  7 Virtual (BPF) control signal u x (t) and u y (t) (blue) and actual control signal v r (t) and ω r (t) (red). The dashed lines represent the physical and the virtual constraints Furthermore, it has been applied successfully to the control of a nonlinear system for which a linearizing/decoupling controller can be devised.
Last, but not by relevance, the proposed approach generates a control signal of the BPF type, that needs, for its generation in a practical implementation, a much simpler controller, making it appealing in the industrial applications.
The current research aims at two goals: (i) to extend the proposed approach to MIMO systems, and (ii) to propose a suitable optimization algorithm to replace the naive bisection search.
Acknowledgements The authors sincerely thank the referees for the constructive comments to the early version of this paper. They also wish to thank their former colleague and mentor professor Luciano Carotenuto of Università della Calabria, Italy, for the helpful discussions.
Author Contributions All authors contributed equally to this work. All authors have read and agreed to the proofs of the manuscript.
Funding Open access funding provided by Universitá della Calabria within the CRUI-CARE Agreement.

Code Availability
The Matlab code for the numerical experiments is available upon request.

Conflicts of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copy-right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.
in particular, the second and the third entries are: In our case we have {h p = h( pT ), p = 1, . . . , m}, where the function h is given by Eq. (4), so Toeplitz's matrix is: We underline that Proposition 1 of Sect. 3 holds for any function h of class C q , q ≥ 0, on the interval [0; mT ].

Proof of Proposition 1
As {h p = h( pT )}, the direct computation of the first b p 's gives , ... .

Now two cases arise.
Case 1: h(0) = 0. In this case, for sufficiently small T we can write h( pT ) ≈ h(0) + p h (1) (0) T , p = 1, . . . , m; as a  consequence, entries b 1 , b 2 If T is small enough that h (1) (0) T can be neglected with respect to h(0), we have and it is easily shown by induction that all the elements b p vanish for p ≥ 3. As a consequence, matrix B T is a bidiagonal matrix that does not depend on T , and its infinity norm is approximated by which proves the thesis for Case 1 with the position η(0) = 2/|h(0)|. Case 2: h(0) = 0 and q is the order of the first derivative that does not vanish at the origin. In this case, for sufficiently small T we can write h( pT ) ≈ h (q) (0) (pT ) q q! , p = 1, . . . , m, and h( pT ) h(T ) ≈ h (q) (0) (pT ) q q! q! h (q) (0) T q = p q , thus, the approximation of the first b p 's yields The structure of those elements lets us conjecture that where γ : N 2 → Z. This is true for p = 1 with γ (1, q) = −1; using the induction principle, we suppose that this is true up to p − 1, and then we compute the generic b p by Eq. (17), obtaining: thus our conjecture is proved by the induction property on defining recursively γ (1, q) = −1, It can also be shown by induction that if q = 1 we get b p ≈ 0 for all p ≥ 4, and if q = 2 we obtain b p+1 ≈ −b p for all p ≥ 4.