Model predictive control for offset-free reference tracking of fractional order systems

TL;DR

This paper introduces an offset-free model predictive control method for fractional-order systems, utilizing a finite-dimensional approximation and disturbance observer to achieve accurate reference tracking and disturbance rejection.

Contribution

It presents a novel offset-free MPC scheme for fractional-order systems using the Grünwald-Letnikov derivative with a disturbance observer for steady-state error correction.

Findings

Effective disturbance reconstruction at steady state

Improved reference tracking in fractional pharmacokinetic models

Demonstrated success in controlling anti-arrhythmic medication

Abstract

In this paper an offset-free model predictive control scheme is presented for fractional-order systems using the Gr\"unwald-Letnikov derivative. The infinite-history fractional-order system is approximated by a finite-dimensional state-space system and the modeling error is cast as a bounded disturbance term. Using a state observer, it is shown that the unknown disturbance at steady state can be reconstructed and modeling errors and other persistent disturbances can be attenuated. The effectiveness of the proposed controller-observer ensemble is demonstrated in the optimal administration of an anti-arrhythmic medicine with fractional-order pharmacokinetics.

Figures

Tables

Table 1. Table 1: Values of the performance indicator J 𝐽 J for perturbed values of pharmacokinetic parameters for N = 60 𝑁 60 N=60 and ν = 25 𝜈 25 \nu=25 .

	$-10\%$	Nominal	$+10\%$
$k_{10}$	9.9628	11.6345	12.7328
$k_{12}$	11.1457		12.0466
$k_{21}$	12.0169		11.2554
$\alpha$	12.4379		10.7575

Table 2. Table 2: Values of the performance indicator J 𝐽 J for different values of the memory length for N = 60 𝑁 60 N=60 .

$\nu$	$J$
5	11.4950
15	11.5998
25	11.6345
35	11.5931

Equations

$${\Delta }_h^{\alpha }f(t)=\sum _{j=0}^{\infty }(-1{)}^j(\genfrac{}{}{0ex}{}{\alpha }{j})f(t-jh)$$

$${\mathrm{D}}^{\alpha }f(t)=\underset{h\to 0}{\lim }\frac{{\Delta }_h^{\alpha }f(t)}{h^{\alpha }}$$

$$\sum _{i=0}^l{A}_i{\mathrm{D}}^{{\alpha }_i}x(t)=\sum _{i=0}^r{B}_i{\mathrm{D}}^{{\beta }_i}u(t)$$

$$\sum _{i=0}^l{\bar{A}}_i{\Delta }_h^{{\alpha }_i}{x}_{k+1}=\sum _{i=0}^r{\bar{B}}_i{\Delta }_h^{{\beta }_i}{u}_k$$

$${\Delta }_{h,\nu }^{\alpha }{x}_k=\sum _{j=0}^{\nu }{c}_j^{\alpha }{x}_{k-j},$$

$$\sum _{i=0}^l{\bar{A}}_i{\Delta }_{h,\nu }^{{\alpha }_i}{x}_{k+1}=\sum _{i=0}^r{\bar{B}}_i{\Delta }_{h,\nu -1}^{{\beta }_i}{u}_k,$$

$$\sum _{j=0}^{\nu }{\hat{A}}_j{x}_{k-j+1}=\sum _{j=0}^{\nu }{\hat{B}}_j{u}_{k-j}$$

$$x_{k+1}=\sum _{j=0}^{\nu -1}{\tilde{A}}_j{x}_{k-j}+\sum _{j=1}^{\nu }{\tilde{B}}_j{u}_{k-j}+{\tilde{B}}_0{u}_k$$

$${\tilde{x}}_{k+1}=A{\tilde{x}}_k+B{u}_k.$$

$${\Delta }_h^{\alpha }={\Delta }_{h,\nu }^{\alpha }+R_{h,\nu }^{\alpha },$$

$$R_{h,\nu }^{\alpha }(x_k)=h^{-\alpha }\sum _{j=\nu +1}^{\infty }{c}_j^{\alpha }{x}_{k-j}.$$

$$R_{h,\nu }^{\alpha }(x_k)\in \underset{j=\nu +1}{\overset{\infty }{⨁}}{c}_j^{\alpha }\mathcal{X}.$$

$$R_{h,\nu }^{\alpha }(x_k)\in \left(\sum _{j=\nu +1}^{\infty }|c_j^{\alpha }|\right)\mathcal{X}.$$

$${\tilde{x}}_{k+1}=A{\tilde{x}}_k+B{u}_k+G{d}_k$$

$${\mathcal{D}}_{\nu }={\mathcal{D}}_{\nu }^x\oplus {\mathcal{D}}_{\nu }^u$$

$${\mathcal{D}}_{\nu }^x$$

$${\mathcal{D}}_{\nu }^u$$

$$y_k=C{\tilde{x}}_k+C_d{d}_k.$$

$${\xi }_{k+1}$$

$$y_k$$

$$\bar{A}=\left[\begin{array}{cc}A& G\\ 0& I\end{array}\right],\bar{B}=\left[\begin{array}{c}B\\ 0\end{array}\right],\bar{C}=\left[\begin{array}{cc}C& C_d\end{array}\right]$$

$$\left[\begin{array}{cc}A-I& G\\ C& C_d\end{array}\right],$$

$${\hat{\xi }}_{k+1}$$

$$e_k$$

$${\mathsf{P}}_{\mathrm{mpc}}({\hat{\tilde{x}}}_k,{\hat{d}}_k):\underset{\text{u}=\{u_{k+j\mid k}{\}}_{j=0}^{N-1}}{\mathrm{minimize}}{V}_N(\text{u};{\bar{x}}_k,{\bar{u}}_k)$$

$${\tilde{x}}_{k+j+1\mid k}$$

$$d_{k+j+1\mid k}$$

$$F_x{\tilde{x}}_{k+j\mid k}$$

$${\tilde{x}}_{k\mid k}$$

$$d_{k\mid k}$$

$$V_N(\mathbf{u};{\bar{x}}_k,{\bar{u}}_k)$$

$$+\sum _{j=0}^{N-1}\left(\Vert {\tilde{x}}_{k+j\mid k}-{\bar{x}}_k{\Vert }_Q^2+\Vert u_{k+j\mid k}-{\bar{u}}_k{\Vert }_R^2\right)$$

$$\left[\begin{array}{cc}A-I& B\\ C& 0\end{array}\right]\left[\begin{array}{c}{\bar{x}}_k\\ {\bar{u}}_k\end{array}\right]=\left[\begin{array}{c}-G{\hat{d}}_k\\ r_k-C_d{\hat{d}}_k\end{array}\right].$$

$$\left[\begin{array}{cc}A-I& B\\ C& 0\end{array}\right]W=\left[\begin{array}{cc}-G& 0\\ -C_d& I\end{array}\right].$$

Full text

††journal: Control Engineering Practice

1 Introduction

1.1 Background

Fractional calculus is a mathematical extension of the classic calculus of integer-order derivatives and integrals. In fractional calculus, derivatives and integrals are extended to non-integer orders which possess fascinating properties. One of the most remarkable properties of fractional-order derivatives is that they are nonlocal operators, that is, unlike their integer-order counterparts, they cannot be evaluated at a certain point solely by knowing how the function behaves in a neighborhood of this point; instead, the whole history of the function needs to be taken into account [1].

Fractional dynamics have been used to model phenomena exhibiting hereditary properties and long or infinite memory transients. Such phenomena include semi-infinite transmission lines with losses [2], viscoelastic polymers [3], magnetic core coils [4], ultra capacitors [5], anomalous diffusion in semi-infinite transmission bodies [6] and several bio-medical applications [7, 8, 9]. Fractional systems find also several applications in physics [3]. Podlubny et al. [10] offer a thorough overview of the wealth of available toolboxes and software that allow the simulation and controller design for fractional-order systems. Fractional-order systems and controllers have made their appearance in industrial applications by the extension of the classical PID controller to fractional-order PIλDμ ones [11, 12, 13, 14]. In [15] fractional order control is applied for the attenuation of vibrations in flexible structures.

During the last few years, a number of works appeared in the literature on the development of Model Predictive Control (MPC) methodologies for fractional order systems. MPC has gained great popularity in industry and academia due to its inherent capability to take into account state and input constraints, handle complex system dynamics and be resilient to external disturbances [16]. In [17, 18] and [19], MPC formulations were presented, based on simple input-output fractional order models and using integer-order approximations of the transfer function of the system. In [19], the proposed fractional order MPC was demonstrated on the low-speed control of gasoline-propelled cars. In [20], both input-output and state space fractional order models were considered as predictive models in MPC. In [21] the use of fractional order Takagi-Sugeno fuzzy models was proposed in the synthesis of fractional MPC.

A common limitation of all aforementioned works is that although they recognize the importance of MPC in handling input and output/state constraints, they do not take explicitly into account those constraints in the proposed MPC formulations. In [22] input constraints were considered in the formulation of the MPC optimization problem, however, without constraints on the state variables and without theoretical stability guarantees.

In practical applications, integer-order derivatives are often used to approximate fractional-order systems [23]. Unfortunately these methods come without guarantees of stability and satisfaction of constraints. In particular, when these approximations are based on frequency-domain procedures, no stability guarantees can be derived and this is a severe shortcoming in safety-critical applications such as drug administration. It should be noted again that the behavior of fractional-order systems depends on the whole history of their trajectories, therefore it is very difficult to provide cogent evidence based on simulations alone without theoretical backup.

An alternative approach has been proposed by Guermah et al. [24] where fractional-order systems are modeled as infinite-dimensional state space systems leading to theoretically interesting stability conditions which, nonetheless, are not tractable and cannot be used for controller design.

In our previous work, a controller design approach based on the Grünwald-Letnikov scheme [25, 26, 9] was proposed. A finite-dimensional approximation was introduced to arrive at a linear time-invariant system and cast the discrepancy between the real and the approximate system as an additive bounded uncertainty term. A worst-case MPC formulation was presented which leads to asymptotically stable behavior towards the origin in presence of state and input constraints, even when an approximate finite-history model is employed, unlike alternative approaches [27, 24].

1.2 Contributions

In this paper an MPC formulation is proposed that achieves offset-free reference tracking for fractional-order systems taking into account input/state constraints. Modeling error is cast as a disturbance term with which the state-space model of the nominal dynamical system is augmented. A state observer is then used to simultaneously estimate the system state and the disturbance using a simple disturbance model. As a result, the closed-loop system can reject disturbances associated with the aforementioned finite-memory approximation, but also other modeling errors due to inexact knowledge of the system parameters, while guaranteeing constraint satisfaction.

Unlike the controller design approaches we discussed in Section 1.1, which use frequency-domain approximations (i.e., integer-order approximations of the transfer function), the approach we propose in this paper uses a time-domain approximation based on the Grünwald-Letnikov derivative.

The proposed MPC strategy is demonstrated in a case study emerging from pharmacokinetics and pharmacodynamics, where fractional-order systems are becoming increasingly popular over the last years. The work of Kytariolos et al. [28] introduced fractional-order dynamics in pharmacokinetics, highlighting why the classical in-vitro-in-vivo correlations theory fails. Certain nonlinearities, anomalous diffusion, deep tissue trapping, diffusion across fractal manifolds such as systems of capillaries, synergistic and competitive actions are cases that can hardly be modeled by integer-order systems [29]. Fractional-order pharmacokinetic dynamics can be cast as physiologically based pharmacokinetic (PBPK) or compartmental models by properly re-writing the mass balance equations using fractional-order derivatives in a way that mass balances are not violated [30]. The controlled drug administration for drugs with fractional dynamics is a key enabler of an effective and realistic therapy and a valuable tool for the clinical practice, especially in presence of constraints [31].

This paper is organized as follows: In Section 2 we describe fractional-order dynamical systems in terms of the Grünwald-Letnikov derivative and derive control-oriented approximations of bounded error. In Section 3 we propose an MPC scheme for offset-free control using a state observer for an augmented system which is able to attenuate modeling errors. Lastly, In Section 4 we present such an offset-free MPC for the control of amiodarone administration to patients and show that the controlled system is resilient to inexact knowledge of the pharmacokinetic parameters of the patients (which are, typically, not known).

1.3 Notation

Hereafter, $\mathbb{R}$ and $\mathbb{N}$ denote the sets of real and nonnegative integers respectively. We denote by $\mathbb{N}_{[k_{1},k_{2}]}$ the set of all integers in the closed interval $[k_{1},k_{2}]$. The set of real $n$-dimensional vectors is denoted by $\mathbb{R}^{n}$ and the set of $m$-by-$n$ matrices by $\mathbb{R}^{m\times n}$. All sets are denoted by calligraphic uppercase letters and all matrices are denoted by uppercase letters. Vectors and scalars are denoted by lowercase letters. The transpose of a matrix $A$ is denoted by $A^{\prime}$.

2 Fractional-Order Systems

2.1 Discrete-time fractional operators

In this section a fractional-order differential operator, the Grünwald-Letnikov derivative is introduced. Let $f:\mathbb{R}\to\mathbb{R}^{n}$ be a bounded function. Let us first introduce the Grünwald-Letnikov difference of $f$ at $t$ of order $\alpha>0$ and step size $h>0$, which is defined as

[TABLE]

where $\binom{\alpha}{0}=1$ and for $j\in\mathbb{N},j>0$, $\tbinom{\alpha}{j}=\prod_{i=0}^{j-1}\frac{\alpha-i}{i+1}$. Furthermore, let us define $c_{j}^{\alpha}=(-1)^{j}\tbinom{\alpha}{j}$ and note that $|c_{j}^{\alpha}|\leq\alpha^{j}/j!$ for all $j\in\mathbb{N}$, therefore, the sequence $c_{j}^{\alpha}$ is absolutely summable and $\Delta_{h}^{\alpha}$ is well defined. It is now clear that in order to estimate $\Delta_{h}^{\alpha}f(t)$ for non-integer orders $\alpha$ the whole history of $f(t)$ is needed.

The Grünwald-Letnikov operator leads to the definition of the Grünwald-Letnikov fractional derivative of order $\alpha$ as

[TABLE]

provided that the limit exists. It can be verified that $\mathrm{D}^{k}$, for $k\in\mathbb{N}$, are the ordinary integer-order derivatives with respect to $t$ and, by convention, $\mathrm{D}^{0}f(t)=f(t)$.

Using the definition above, it is easy to describe fractional-order dynamical systems with state $x:\mathbb{R}\rightarrow\mathbb{R}^{n}$ and input $u:\mathbb{R}\rightarrow\mathbb{R}^{m}$ as

[TABLE]

where $l,r\in\mathbb{N}$, $A_{i}$ and $B_{i}$ are matrices of appropriate dimensions and all powers $\alpha_{i}$ and $\beta_{i}$ are nonnegative.

For the discretization of a fractional system, an Euler-type method is used to approximate $\mathrm{D}^{\alpha}$ in (3), for a fixed time step size $h$, using $h^{-\alpha}\Delta_{h}^{\alpha}$. Using the forward operator for the derivatives of the states, and the backward operator for the input variables, the discretization of Equation (3) becomes

[TABLE]

For convenience in (4) it is: $x_{k}=x(kh),u_{k}=u(kh)$ for $k\in\mathbb{Z}$, and $\bar{A}_{i}=h^{-\alpha_{i}}A_{i},\bar{B}_{i}=h^{-\beta_{i}}B_{i}$.

As it becomes obvious from the equations above, in order to fully describe a fractional-order system, infinite-dimensional operators of the form $\Delta_{h}^{\alpha}$ should be used. Therefore, it is practically impossible to simulate such systems or design feedback controllers. In what follows, Equation (4) will be approximated using the methodology described in [25] where a finite-dimensional approximation is proposed and the approximation error is treated as a bounded additive disturbance.

2.2 Bounded-error Finite-dimension Approximation

As infinite-memory systems pose severe limitations regarding their simulation and controller design, the truncated Grünwald-Letnikov difference operator of length $\nu$ is introduced

[TABLE]

and the approximation of (3) for $\nu\geq 1$ becomes

[TABLE]

which is equivalently written as

[TABLE]

where $\hat{A}_{j}=\sum_{i=0}^{l}\bar{A}_{i}c_{j}^{\alpha_{i}}$ and $\hat{B}_{j}=\sum_{i=0}^{r}\bar{B}_{i}c_{j}^{\beta_{i}}$ for $j\in\mathbb{N}_{[0,\nu]}$. Assuming that matrix $\hat{A}_{0}$ is nonsingular, $\tilde{A}_{j}=-\hat{A}_{0}^{-1}\hat{A}_{j+1}$ and $\tilde{B}_{j}=\hat{A}_{0}^{-1}\hat{B}_{j}$ are defined and the last Equation (7) becomes

[TABLE]

Considering as the state variable $\tilde{x}_{k}=(x_{k},x_{k-1},...,x_{k-\nu+1},u_{k-1},...,u_{k-\nu})^{\prime}$, the system can be written in the compact form of a finite-dimensional linear time-invariant (LTI) system as

[TABLE]

In the next section, this form will be used to formulate a model predictive control problem, assuming that the pair $(A,B)$ is stabilizable.

The fractional-order difference operator can be written as the sum of the truncated difference operator and a residual term containing all those past states that were not taken into account, that is

[TABLE]

where $R_{h,\nu}^{\alpha}$ is the operator

[TABLE]

Assume that up to time $k$ the state has been constrained in a compact convex set $\mathcal{X}\subseteq\mathbb{R}^{n}$ containing the origin in its interior, that is $x_{k-j}\in\mathcal{X}$ for all $j\in\mathbb{N}$. Let us denote the Minkowski sum of two sets $\mathcal{C},\mathcal{D}\subseteq\mathbb{R}^{n}$ as the set $\mathcal{C}\oplus\mathcal{D}\mathrel{:}=\{c+d:c\in\mathcal{C},d\in\mathcal{D}\}\subseteq\mathbb{R}^{n}$. For a collection of sets $\mathcal{C}_{1},\mathcal{C}_{2},\ldots$, their Minkowski sum $\mathcal{C}_{1}\oplus\mathcal{C}_{2}\oplus\cdots$ is denoted as $\bigoplus_{k=1}^{\infty}\mathcal{C}_{k}$.

By virtue of (11),

[TABLE]

In the special case where $\mathcal{X}$ is also a balanced set, that is $-x\in\mathcal{X}$ for all $x\in\mathcal{X}$, (12) reduces to

[TABLE]

In this case, it is evident that $\nu$ directly controls the size of the right hand side of (13). In particular, $({}\sum_{j=\nu_{1}+1}^{\infty}|c_{j}^{\alpha}|{})\mathcal{X}\supset({}\sum_{j=\nu_{2}+1}^{\infty}|c_{j}^{\alpha}|{})\mathcal{X}$ for all $\nu_{1}<\nu_{2}$, that is, the larger the value of $\nu$ is, the lower the worst-case approximation error will be.

The approximation error can be integrated into the state space form as an additive disturbance term, so (8) becomes

[TABLE]

with $G=\begin{bmatrix}I&0&\cdots&0\\ \end{bmatrix}^{\prime}$. The disturbance term $d_{k}$ is introduced to cast the approximation error which is in fact a bounded error term as we are about to show. Assuming that the state is constrained, as discussed above, in a set $\mathcal{X}$ and, similarly, the input is constrained in a convex compact set $\mathcal{U}$ containing the origin in its interior, following [25], $d_{k}$ is bounded in a compact set $\mathcal{D}_{\nu}$ given by

[TABLE]

where sets $\mathcal{D}_{\nu}^{x}$ and $\mathcal{D}_{\nu}^{u}$ are111Simpler formulas for $\mathcal{D}_{\nu}^{x}$ and $\mathcal{D}_{\nu}^{u}$ may be derived assuming that $\mathcal{X}$ and $\mathcal{U}$ are balanced sets, i.e., they are symmetric about the origin — or equivalently $-x\in\mathcal{X}$ whenever $x\in\mathcal{X}$ and $-u\in\mathcal{U}$ whenever $u\in\mathcal{U}$. The reader is referred to [25] for further information.

[TABLE]

In the special case where $\mathcal{X}$ and $\mathcal{U}$ are balanced sets, these expressions can be further simplified as explained above.

The assumption that $\mathcal{X}$ and $\mathcal{U}$ contain the origin in their interiors is standard in regulation problems, where the objective is to steer the state of the system to the origin. In case of reference tracking, it is desirable that the system output $y_{k}$, which is typically a linear combination of states $y_{k}=Cx_{k}$, converges to a certain set-point $r$. Assume that the set-point is attainable by a feasible state $\bar{x}\in\mathcal{X}$ so that $C\bar{x}=r$. It then suffices that $x_{k}$ converges to $\bar{x}$ with $x_{k}\in\mathcal{X}$ with a feasible sequence of control actions $u_{k}\in\mathcal{U}$ which converges to a feasible point $\bar{u}\in\mathcal{U}$. Then, the deviation variables $\delta x_{k}=x_{k}-\bar{x}$ and $\delta u_{k}=u_{k}-\bar{u}$ may be defined. The above convergence requirements become $\delta x_{k}\to 0$ and $\delta u_{k}\to 0$ with constraints $\delta x_{k}\in\{-\bar{x}\}\oplus\mathcal{X}$ and $\delta u_{k}\in\{-\bar{u}\}\oplus\mathcal{U}$. In order to be able to derive bounded-error approximations of the fractional order system with input $\delta u_{k}$ and state $\delta x_{k}$ it needs to be assumed that $\bar{x}$ lies in the interior of $\mathcal{X}$ and $\bar{u}$ lies in the interior of $\mathcal{U}$.

Larger memory lengths lead to better approximations since $R_{h,\nu}^{\alpha}(x)\to 0$ as $\nu\to\infty$, that is, arbitrarily low approximation errors can be obtained by choosing an appropriately large value of $\nu$. That said, on one hand it is desirable to opt for a large value of $\nu$ to obtain better approximations. On the other hand, too large values of $\nu$ should be avoided as they will incur a larger computation cost especially regarding the solution of the MPC problem in real time (see Section 3.2).

3 Offset-free reference tracking

3.1 Observer design

In this section, it is assumed that the output of (14) can be measured, which is given by

[TABLE]

A state observer will be designed to simultaneously reconstruct the state $\tilde{x}_{k}$ and persistent disturbances $d_{k}$ using the simple disturbance model $d_{k+1}=d_{k}$. In [32] it is shown that even if the disturbance model is not accurate, the offset-free regulation properties of the closed-loop system and the closed-loop performance will not be altered.

We introduce an augmented state observer which is a dynamical system which using input-output measurements, will produce estimates of the augmented state variable \xi_{k}=\big{[}\tilde{x}_{k}^{\prime},d_{k}^{\prime}\big{]}^{\prime} whose dynamics is

[TABLE]

with

[TABLE]

where $I$ is the identity matrix. According to [33, Prop. 1], system (18) is observable if and only if $(C,A)$ is observable and

[TABLE]

is full column rank. For these conditions to be satisfied, we need to choose the dimension of $d_{k}$ to be no larger than the number of outputs. We shall hereafter assume that $d_{k}$ has the same dimension as $y_{k}$. Additionally, matrix $C_{d}$ reflects the effect that the persistent disturbance $d_{k}$ has (directly) on the system output.

The observer is a dynamical system that produces estimates \hat{\xi}_{k}=\big{[}\hat{\tilde{x}}_{k}^{\prime},\hat{d}_{k}^{\prime}\big{]}^{\prime}, that is, the observer produces simultaneously estimates of the system state and the disturbance. Here we use a linear observer of the form

[TABLE]

where $e_{k}$ is the state estimation error.

By properly choosing the gain matrix $L$, either by pole placement methods or with LQG methods, it can be guaranteed that the estimation error converges to [math] for any initial estimate $\hat{\xi}_{0}$. The estimates of the augmented state are then provided to the tracking MPC which is described in the following section.

3.2 Offset-free model predictive control

In MPC, the control action is computed at every time instant by minimizing an index which quantifies the performance of the system along a finite prediction horizon. This performance index is used to choose an optimal sequence of control actions among the set of such admissible sequences, while corresponding state sequences are predicted using a system model. The first element of the optimal sequence is then applied to the system.

Maeder et al. [33] proposed an offset-free MPC problem formulation without integral action where the MPC problem makes use of the augmented model (18). The main concept of this approach is that the observer eventually reconstructs any persistent disturbances and this information is used to steer the state to the desired set-point ${r_{k}}$ without offset. The objective of the closed-loop system is for the tracking error $y_{k}-r_{k}$ to converge to [math] whenever $r_{k}$ is a convergent sequence.

At every time instant $k$ the following finite-horizon optimal control problem is solved:

[TABLE]

subject to the constraints

[TABLE]

with cost function

[TABLE]

where the notation $\|x\|_{F}^{2}=x^{\prime}Fx$ is used. The optimization problem (22) is solved over sequences of future control actions $\textbf{u}=\{u_{k+j\mid k}\}_{j=0}^{N-1}$ of length $N$ and the problem is solved using information that is available at time $k$ which is provided by the augmented state observer. The first element of this sequence, $u_{k\mid k}^{\star}$, is then applied to the system in a receding horizon fashion [16].

The optimization problem $\mathsf{P}_{\mathrm{mpc}}(\hat{x}_{k},\hat{d}_{k})$ is solved with respect to the dynamics of the augmented system described by (23a) and (23b), state and input constraints (23c) and the initial conditions (23d) and (23e). The constraints (23c), where $F_{x}$ and $F_{u}$ are matrices of appropriate dimensions and $\leq$ is meant element-wise, encompasses simple bounds of the form $\tilde{x}_{\min}\leq\tilde{x}_{k}\leq\tilde{x}_{\max}$ and $u_{\min}\leq u_{k}\leq u_{\max}$ or more involved polyhedral joint state-input constraints.

The terms $\bar{x}_{k}$ and $\bar{u}_{k}$ in (22) are calculated at runtime by solving the linear system

[TABLE]

Equation (25) is used at every time instant to determine $\bar{x}_{k}$ and $\bar{u}_{k}$ given the desired set-point $r_{k}$ and the disturbance estimate $\hat{d}_{k}$. In fact, it suffices to find a matrix $W$ so that

[TABLE]

Then, $\bar{x}_{k}$ and $\bar{u}_{k}$ can be obtained by

[TABLE]

As the linear observer produces a sequence of disturbances $\hat{d}_{k}$ which converges to a steady state value $\hat{d}_{\infty}=d_{\infty}$, conditions (25) enforce that the tuple ($\bar{x}_{\infty},d_{\infty},\bar{u}_{\infty},r_{\infty})$ is an equilibrium point of (18), i.e., as $r_{k}\to r_{\infty}$, $x_{k}\to\bar{x}_{\infty}$, $u_{k}\to\bar{u}_{\infty}$, we have $y_{k}\to r_{\infty}$.

Finally, the positive semidefinite matrix $Q$ and the positive definite matrix $R$ are the tuning knobs of MPC and are typically chosen to be diagonal matrices; larger values of $Q$ lead to a faster convergence of the system states to their equilibria, whereas larger values of $R$ lead to a smoother actuation and slower convergence towards the steady state. Matrix $P$, which defines the terminal cost function in (3.2), is taken to be the unique solution of the Riccati equation

[TABLE]

By choosing $P$ to be the solution of (28), the value function of (22) becomes a Lyapunov function for the closed-loop system and stability is guaranteed [33, 16].

The optimization problem $\mathsf{P}_{\mathrm{mpc}}$ is a convex quadratic optimization problem with polyhedral constraints which can be solved very efficiently in practice. In fact, there exist algorithms which exhibit very fast convergence and scale linearly with the prediction horizon [34, 35].

4 Application: Controlled drug administration

4.1 Problem statement: objectives and constraints

Fractional-order dynamics are of great interest in applications of pharmacokinetics such as controlled drug administration [30, 36, 31]. In this section the reference tracking methodology presented in Section 3 is applied for the controlled intravenous administration of a drug which exhibits fractional-order dynamics.

Amiodarone (CAS registration number: 1951-25-3) is an anti-arrhythmic medication used for the treatment of ventricular tachycardia, shock-refractory ventricular and atrial fibrillation, which has been successfully modeled by a fractional compartmental system in [30]. The pharmacokinetic model considers two compartments: the central compartment which corresponds to the blood stream where the drug is introduced with rate $u$ and a peripheral compartment which corresponds to the tissues where the blood is distributed. The distribution from the blood stream to the tissues is assumed to follow a simple first-order dynamics with constant $k_{12}$. Amiodarone is excreted from the central compartment following first-order excretion kinetics with rate constant $k_{10}$. A mass flow from the tissues back to the central compartment takes place via anomalous diffusion and follows fractional-order kinetics of order $1-\alpha$ with rate $k_{21}$. This structure is shown in Figure 1 along with the observer and the offset-free controller.

The mass balance equations, taking into account the fractional-order kinetics, lead to the dynamical system

[TABLE]

where $A_{1}$ and $A_{2}$ are the amounts of amiodarone (in $\mathrm{ng}$) in the central and peripheral compartments respectively and $u$ is the intravenous administration rate into the central compartment (in $\nicefrac{\mathrm{ng}}{\mathrm{day}}$). The (nominal) values of the parameters of the system are $\alpha=0.587$, $k_{10}=1.4913\,\mathrm{day^{-1}}$, $k_{12}=2.9522\,\mathrm{day^{-1}}$ and $k_{21}=0.4854\,\mathrm{day^{-\alpha}}$. It is considered that amiodarone is administered to the patient intravenously and continuously. Only measurements of $A_{1}$ are available in real time.

As shown in Figure 1, a state observer receives measurements of the system output, that is $y=A_{1}$ and input variables and produces estimates of the full state $\hat{\tilde{x}}$ and disturbance $\hat{d}$ which are provided to the MPC. The treating physician prescribes the desired set-point $r$ for $A_{1}$ to the MPC controller which decides the administration rate $u$ at every time instant.

The continuous-time fractional-order system is simulated using the Oustaloup filter — an integer-order approximation of the fractional-order differentiators $s^{a}$ in a frequency range $[\omega_{L},\omega_{H}]$ via a logarithmically spaced sampling of $N_{f}$ points [37, 38]. Using the Oustaloup filter, $s^{a}$ is approximated by a transfer function with numerator and denominator degrees equal to $N_{f}$.

Oustaloup’s filter is chosen because, as shown in [39], it is very accurate for this dynamical system. Following [39], the values $\omega_{L}=10^{-2}\,\mathrm{day^{-1}}$, $\omega_{H}=10^{3}\,\mathrm{day^{-1}}$ and $N_{f}=8$ are chosen which are known to lead to a highly accurate approximation. Other frequency-domain approximation methods are available in the literature such as the Padé and the Matsuda-Fujii which are studied in [39] in terms of accuracy and seem to lead to solutions of comparable quality. We should underline, however, that frequency domain methods do not come with bounded error guarantees, that is, we cannot know the maximum approximation error in advance. This motivates the use of the Grünwald-Letnikov approximation which is the most suitable approach. Frequency-domain approximations lead to approximation errors which — although in practice are small — are not accompanied by known bounds. On the other hand, the truncated Grünwald-Letnikov approximation, as less parsimonious as it might be, leads to a bounded-error approximation in the time domain which is necessary for the design of MPC controllers with theoretical stability and constraint satisfaction guarantees. Other methods, such as the numerical inverse Laplace transform, lead to highly accurate solutions, but these are not suitable for controller design.

In practice people have used other approximations to design MPC controllers (than the Grünwald-Letnikov discrete-time derivative) [20, 40, 22], however, without any guarantees of constraint satisfaction which is particularly important in fail-sensitive applications (e.g., medical).

4.2 Control-oriented modeling

The fractional-order pharmacokinetic system (29) can be written as

[TABLE]

with $\beta=1-\alpha$ and $k_{0}=k_{12}+k_{10}$. The constants

[TABLE]

are defined and the variable

[TABLE]

is introduced. Now (30) can be rewritten in a compact form as

[TABLE]

The above system can now be discretized with sampling time $h$ as discussed in Section 2.2 to yield

[TABLE]

where $B_{d}=[\begin{smallmatrix}1\\ 0\end{smallmatrix}]$ and $d_{k}\in\mathbb{R}$ is a generic disturbance term which is used to encompass the effect that the following factors have on $A_{1}$: (i) the approximation error due to the use of the truncated operator $\Delta_{h,\nu}^{\beta}$, (ii) modeling errors due to inexact knowledge of the actual pharmacokinetic parameters, (iii) modeling errors due to pharmacokinetic dynamics not captured by the model described in the the previous section. We shall now use Equation (34) to derive a simplified linear time-invariant model with a scalar disturbance term which will be used to design the MPC controller and the observer as discussed above.

Let us define the state variable $\tilde{x}_{k}=(x_{k},x_{k-1},\ldots,x_{k-\nu})\in\mathbb{R}^{2(\nu+1)}$. By defining $\Lambda=Mh+I+\Theta h^{\beta}$ the approximate system in a matrix form is derived:

[TABLE]

therefore, $C=\begin{bmatrix}1&0&\ldots&0\end{bmatrix}$ and $C_{d}=0$. One may also verify that the observability assumptions stated in Section 3.1 are satisfied for system (35).

In what follows the initial condition is $A_{1}(t)=A_{2}(t)=0$ for all $t\leq 0$ is assumed, that is $\tilde{x}_{0}=0$.

4.3 Controlled drug administration

Hereafter a discretization with $h=0.1\,\mathrm{day}$ and a memory length $\nu=25$ which corresponds to $2.5$ days are used. In the closed-loop simulation, we assume that the treating physician sets the reference to the system output ($A_{1}$) equal to $0.5\,\mathrm{ng}$ for the first $80$ days and increases it to $1.0\,\mathrm{ng}$ for the next $70$ days. The administration rate is constrained in the interval $0\leq u_{k}\leq u_{\max}=2\,\nicefrac{\mathrm{ng}}{\mathrm{day}}$, while the system output should not exceed the upper bound of $1.03\,\mathrm{ng}$ . The prediction horizon for the controller is set to $N=60$. The weight matrices of the MPC controller are chosen to be $Q=0.25\cdot I$ and $R=5$.

The quantity of amiodarone in the central compartment, $A_{1}$, is shown in Figure 2 (left) where it can be seen that it tracks the prescribed set-point while it does not exceed the constraint of $1.03\,\mathrm{ng}$. The input produced by the controller as the solution of the optimization problem $\mathsf{P}_{\mathrm{mpc}}$ is shown in Figure 2 (right) and it may be observed that the administration rate does not exceed the maximum allowed value of $2.0\,\nicefrac{\mathrm{ng}}{\mathrm{day}}$. In Figure 2 (middle) we show the amount of drug in the peripheral/tissues compartment.

4.4 Sensitivity analysis

In order to assess the sensitivity of the closed-loop system to the exact knowledge of the pharmacokinetic parameters $k_{12}$, $k_{21}$, $k_{10}$ and $\alpha$ we performed closed-loop simulations against perturbed models using the MPC designed for the exact model. In particular we assumed a perturbation of $+10\%$ or $-10\%$ on each one of the pharmacokinetic parameters.

In each closed-loop simulation, we assumed the same scenario as in the nominal case, i.e., the treating physician selects a set-point equal to $0.5\,\mathrm{ng}$ for the first $80$ days and increases it to $1.0\,\mathrm{ng}$ for the next $70$ days. We also assumed the same upper bounds on the administration rate and the quantity of amiodarone in the central compartment, $A_{1}$. In order to evaluate the performance of the drug administration course, we introduced the index

[TABLE]

where $N_{u}$ is the number of simulation points (here $N_{u}=3000$). This performance indicator quantifies and encodes the average tracking error and the average administration rate, both of which are desirable to be low.

Figures  3 and  4 present the input-output responses and the disturbance estimates $\hat{d}_{k}$ and Table  1 shows the value of the performance index for the nominal case and all eight perturbations. It seems that the closed-loop system is more sensitive to errors in the clearance rate constant $k_{10}$, however its effect is rather limited. Evidently, in all simulation results, the impact of parametric uncertainty is not significant and the offset-free control scheme produces successful control actions in tracking the prescribed set-point changes. Overall, despite the fact that the exact values of $k_{10}$, $k_{12}$ and $k_{21}$ are not available to the controller, the prescribed set-point on $A_{1}$ is reached and the constraints on $A_{1}$ and $u$ are satisfied at all time instants.

In the following simulations, we investigated the effect of the value of the memory length $\nu$ on the performance of the system. Again, we ran the same scenario with the nominal case concerning the set point changes and the constraints of the input and state parameters using memory lengths equal to $5$, $15$, $25$ and $35$. As we discussed in Section 2.2, larger values of $\nu$ lead to more accurate models. As one may observe in Figure 5 and Table 2, the performance is only slightly improved by increasing the value of a memory length $\nu$ and in general the closed-loop performance is not sensitive to the selection of $\nu$. For $\nu=25$ the response curves are almost indistinguishable from those with $\nu=35$.

Overall, we see that parametric errors can be well attenuated and have little effect on the performance of the closed-loop controlled system. As one may observe in Figures 3, 4 and 5, the disturbance estimates $\hat{d}_{k}$ reconstruct the model-system mismatch and tracking offset is avoided.

4.5 Computational aspects

All simulations were run in MATLAB® (version 2014b) on a Macbook Pro (version 10.12.4, 2.66 GHz Intel Core 2 Duo, 4 GB RAM) running Mac OS Sierra. The MPC problem was formulated as a constrained quadratic program which can be solved very efficiently in real time. The solver used for solving the optimization problem was Mosek [41]. The average runtime to solve the MPC problem was 1.3 s and the maximum runtime was 2.3 s for a memory of $\nu=35$.

5 Conclusions

In this paper an offset-free tracking MPC scheme was proposed for constrained fractional-order systems and it was demonstrated that it is resilient to modeling errors.

The proposed tracking methodology is based on a discrete-time finite-dimension approximation of the infinite-dimensional fractional dynamics and can handle constraints on both state and input variables. The approximation used in this work hinges on the Grünwald-Letnikov fractional-order derivative, it is a time-domain approximation and, unlike frequency-domain approximations, leads to a bounded approximation error with known error bounds. Combined with the design of an observer, the proposed offset-free MPC approach is compatible with the availability of partial state information, which makes it suitable for real control problems. The augmented state observer which is used in this work produces estimates both for the system state and a disturbance term which encompasses modeling errors. Asymptotically, the observer reconstructs both the states and the disturbance term leading to offset-free tracking despite possible modeling errors.

The merits of the proposed MPC scheme were demonstrated on the problem of intravenous administration of amiodarone; a drug which follows fractional-order pharmacokinetics. Using blood stream measurements only, the closed-loop system, which combines an augmented state observer with an MPC controller, leads to zero-offset tracking and is resilient to modeling errors. This characteristic of the proposed methodology deems it suitable for drug administration applications where the actual pharmacokinetic parameters of the patient are unknown. We demonstrated that the sensitivity of the closed-loop system to the errors on the pharmacokinetic parameters is rather low. Additionally, a memory length of $\nu=25$ is adequate for controller design purposes leading to a computationally tractable formulation. The resulting optimization problems are solved very efficiently and do not pose any limitations to the applicability of the controller.

The theory and the results presented in this paper demonstrate that the proposed MPC scheme inherits the key advantages of offset-free MPC to fractional-order systems in terms of satisfying the system constraints and achieving zero steady-state error even in presence of modeling uncertainties. Future research directions will examine extensions of this scheme to formulate robust or stochastic MPC strategies and take into account time delays.