Offset-Free Nonzero Tracking for Nonlinear Impulsive Systems With Application to Biomedical Processes

Impulsive control systems have shown strong potential to represent and address challenging problems, especially in the biomedical field. In recent research, these problems have been tackled with advances in linear impulsive control systems. However, many biomedical applications are better described by nonlinear impulsive models, and therefore, it is necessary to develop analysis tools and control strategies in this context. In the literature, the main properties of nonlinear impulsive control systems have been fully understood, but there is no major development of control strategies. Particularly, there is no substantiation of model predictive control (MPC) strategies maintaining convexity of the optimization problem and closed-loop stability, and there is no control strategy to reduce the offset problem when there are parameter variations, which is a common situation in biological processes. Therefore, the main novelties of this paper are: (i) an MPC formulation extended to nonlinear impulsive systems that addresses non-zero tracking, (ii) the sufficient and necessary conditions to guarantee the stability of the closed-loop system at an equilibrium target, (iii) a comprehensive description of an offset-free MPC to handle low to moderate plant-model mismatches, (iv) the conditions to guarantee offset-free control. Finally, the MPC and offset-free MPC are tested to address the drug administration problem in two biomedical applications: oncolytic virus therapy, to regulate tumor dynamics using doses of oncolytic, and type 1 diabetes treatment, to regulate glycemia using insulin injections. Satisfactory results were obtained in simulation scenarios including parameter variations in nonlinear models that represent the corresponding dynamics.


I. INTRODUCTION
Impulsive Control Systems (ICSs) are a class of systems in which the input action has a short duration in comparison with the sampling time and the dynamics of the system itself. ICSs are characterized by two responses, a free response corresponding to the evolution of the state when the input is zero, and a forced response observed as instantaneous jumps when the control action is applied [1], [2]. ICSs have received increasing attention in the context of industrial process control, biomedical research, and other applications. In biomedical research, ICSs have been useful to represent some therapies for human immunodeficiency virus (HIV) infection [1], [3], type 1 diabetes mellitus [4], [5], and The associate editor coordinating the review of this manuscript and approving it for publication was Wonhee Kim. oncolytic virus therapy [6], showing that this approach offers ways to improve the dosage and schedule of therapies based on control engineering strategies.
Significant progress has been made in formulating ICSs with model predictive control (MPC). MPC is a control strategy that uses a model to predict the system's evolution, and based on that information, MPC computes an optimal control sequence to force the evolution of the system to fulfill some predefined constraints while minimizing a cost function. The first instance of MPC applied to ICSs was developed in [7], where the state is steered to a target zone that does not include the origin. In [8] and [9] a new formulation of zone MPC (ZMPC) was introduced to ensure feasibility in the case of a change of the target by using virtual/intermediary equilibrium variables to steer the state to an equilibrium target. Additionally, some adaptive impulsive VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ MPC formulations were proposed for type 1 diabetes in [5] and [10]. Regarding an MPC formulation for nonlinear ICSs, a first approach was presented in [11] with application to oncolytic virus therapy. This strategy is based on the standard MPC, where the output is steered to a fixed target point. This MPC strategy was not substantiated in terms of feasibility and stability. A common problem in control systems is the inherent mismatch between the prediction model and the actual plant given parameter or model uncertainty. When there is a mismatch, the model used for the predictions does not describe the plant behavior with sufficient accuracy, and thus, the objective is not reached when applying the control action generated by the MPC, but instead, a steady-state offset is produced [12]. To counteract this situation, a more elaborated control strategy is required. In [13]- [16], approaches focused to achieve offset-free tracking for MPC were developed for discrete-time systems. These formulations mainly consist of augmenting the state with a disturbance, and by means of a state estimator obtain information of both the state and the disturbance to correct the prediction model and the target in the MPC formulation. The linear results were extended to discrete-time nonlinear plants in [17] and [18] where the nonlinear model was used as prediction model in the MPC problem, and in [19] where online linearizations of the nonlinear plant were used for prediction and control of wiener systems. In [20], it was demonstrated that the velocity form and state disturbance observer approaches for offset-free MPC are equivalent to particular cases of the disturbance model and observer. A comparison of different offset-free approaches including the disturbance modeling and the analysis of state estimation algorithms coupled with the MPC was presented in [21].
Additional improvements of the offset-free strategy for discrete-time systems can be found in the literature. For instance, in [22] a two-tiered control structure including the linear offset-free MPC strategy was developed to explicitly incorporate the desired closed-loop behavior in the control design. In [23] an adaptive approach for estimating a time-varying disturbance was proposed, and in [24], the authors proposed a semi-infinite program to generate disturbance models for offset-free nonlinear MPC applicable to small to medium scale models. Also, economic MPC formulations including the offset-free feature to account for model uncertainties can be found in [25]- [27]. A different approach for discrete-time systems was developed in [28] using a multiple linear MPC that can handle some amount of plant-model uncertainty, and it was improved with an adaptive integral action to guarantee offset-free tracking while avoiding the use of an observer.
Regarding impulsive systems, the linear offset-free MPC was developed in [29], where observability conditions for the augmented system with a disturbance are established, and a straightforward way to select the disturbance model matrices is introduced based on the properties of ICS. To the best of the authors' knowledge, offset-free MPC formulations for nonlinear ICS have not yet been presented.
In summary, the aspects that are worthy of being developed or further improved are: 1) Biomedical processes are highly nonlinear, and their treatments are better modeled when considering impulsive inputs. Thus, nonlinear impulsive strategies are required. 2) From the literature review, MPC formulations for linear impulsive systems have been proposed to steer the system to nonzero set points or target sets. For nonlinear impulsive systems, there is no substantiation of an MPC strategy. Also, the theoretical substantiation of different nonlinear MPC strategies for discrete-time systems can be found, but these strategies usually use the nonlinear model as a prediction model producing a non-convex optimization problem that is hard to solve. 3) From the literature review, offset-free MPC strategies for linear and nonlinear discrete-time systems have been found. These cover different approaches with and without a state estimator, robust and adaptable formulations. However, only the linear version of the offset-free strategy for impulsive systems has been developed. There are no MPC strategies for nonlinear impulsive systems that address the offset problem caused by a plant-model mismatch. To contribute to the solution of the previous problems, this paper develops two MPC strategies for nonlinear ICSs and provides the conditions to guarantee the stability of the closed-loop system and then for offset-free control. The main merits of this work are listed as follows: 1) An MPC formulation for nonlinear ICS is developed based on the ZMPC that uses (i) an equilibrium artificial reference as a new decision variable, (ii) a dynamic cost function that penalizes the deviation to the artificial equilibrium, and (iii) a final cost function to steer the artificial variable to the actual target set. This strategy has the advantage of allowing the inclusion of a target set instead of a set-point, and it ensures feasibility under any change of the target. 2) For this nonlinear strategy, conditions to stabilize the system at the defined target are presented. This is done by proving that the stability of the continuous trajectory of the ICS follows from the stability of the discrete points of an underlying subsystem using an affine linear approximation at each time state, which guarantees that the optimization problem remains convex. 3) An offset-free MPC strategy is substantiated for nonlinear ICS. This is done by following the general approach of augmenting the system with a disturbance model that represents the plant-model mismatch and ensuring that this new extended impulsive system is observable. Then, the information of the plant-model mismatch is used to correct the prediction model and target of the MPC to achieve zero offset tracking. An additional advantage of the proposed offset-free formulation is the merge of the target calculation problem and the MPC problem in a single one by using the artificial equilibrium variables. 4) To show the effectiveness of both schemes, the MPC and the offset-free MPC (MPC-OF) are applied to the design of candidate protocols for the treatment of two relevant biomedical applications: oncolytic virus therapy (OVT) and Type 1 diabetes (T1D) treatment. OVT is a promising approach to treat cancer using genetic engineering of viruses to destroy tumor cells without infecting or damaging healthy cells [30]- [34]. Improvements in oncolytic viruses delivered promising experimental findings, which motivated clinical trials for the treatment of cancers such as breast cancer [30], [32], [33], and the development of mathematical models to represent the virus and tumor dynamics [6], [30], [35], [36]. However, the antitumor efficacy of oncolytic viruses tends to decrease in animal and human trials due to several challenges such as dose-limiting toxicity, intratumoral viral infection, antiviral immune responses, and the risk of failure to achieve robustness due to the heterogeneous nature of tumors and virus dynamics [30]- [33], [35]. Thus, delivering robust and efficient therapies motivate interdisciplinary research to resolve challenges impairing clinical outcomes. Unlike former protocols in literature, therapies using the offset-free scheme here presented exhibit some robustness to variations in biological rates, and modeling uncertainty such as the impact of stromal cells.
On the other hand, closed-loop control systems for insulin management in subjects with T1D have seen tremendous progress, leading to the development of the artificial pancreas system [37]. The artificial pancreas consists of 3 main components: a continuous glucose monitor sensor, an insulin infusion pump, and a control algorithm to compute the insulin doses required by the patient [5], [38]. Regarding the control strategy, MPC strategies have shown efficiency in simulation and clinical trials [38]- [40], and recent works have also developed formulations under the scope of ICS [4], [10]. Nevertheless, challenges as intra-and inter-subject variations are still to be overcome. In this work, the nonlinear minimal model developed in [41] is used to assess the performance of the MPC strategies here presented under intra-subject variations in insulin sensitivity and insulin degradation rate.
The outline of this paper is as follows: The Materials and methods Section is introduced with the notation and some preliminaries about impulsive systems. Next, the ZMPC and offset-free schemes are developed for nonlinear ICS. In the Results Section the control strategy is applied to propose oncolytic adenovirus dosage for tumor regulation, and insulin dosage for diabetes treatment. Lastly, the conclusions and further work are exposed.

A. NOTATION
The sets of nonnegative integers, reals, n-dimension column vectors, and matrices of dimension n×m are denoted as N, R, R n and R n×m , respectively. I n denotes the identity matrix of dimension n×n. The transpose of a matrix A is represented as A . The convex hull of a collection of sets V k , k = 1, · · · , k is denoted as ch{V 1 , V 2, · · · , V k }. dist V (x) min y∈V y − x , is the distance from a point x to a nonempty closed set V . Given x ∈ R n and a matrix Q ≥ 0 ∈ R n×n , x 2 Q = x Qx.

B. NONLINEAR IMPULSIVE CONTROL SYSTEMS
The class of dynamic systems of interest in this paper consists of a set of nonlinear impulsive control systems (NICSs) with the following representation: where ξ ∈ X ⊆ R n and u ∈ U ⊆ R m denote the system's constrained state and control inputs. These sets are supposed to contain the origin and to be convex, compact sets. Function f : R n → R n denotes the autonomous response of the nonlinear system in absence of a control input. B ∈ R n×m is a matrix that represents the effects of the control input in the instants τ k , k = 1, 2, . . . over the state ξ . z ∈ Y ⊆ R p represents the system's measurable variables, and C ∈ R p×n is the constant matrix of the output equation.
By considering a fixed sampling time T s defined as T s = τ k+1 − τ k in which the impulses are applied, the following discrete counterpart of the NICS can be associated to system (1): where ξ • represents the state ξ (τ − k ) before the impulsive input, and it is also the state available for measurements. ξ • is the state ξ (τ + k ) after applying the impulsive input at times τ k , with k = 0, 1, 2, · · · . f d (·) represents the discretized function of f (·), and φ(ξ, t − τ k ) is the solution of the autonomous part of (1) during the time interval t ∈ [τ k , τ k+1 [, with initial condition ξ . This solution denotes the orbit of ξ defined as o ξ = φ(ξ, t − τ k ), which starts after the impulse at ξ • and finishes at ξ • . It should be noted that each discrete nonlinear difference equation in (2) depends only on its respective state (i.e., only on ξ • or ξ • , respectively) and the input.
Consider the discrete NICS given by (2). A state (ξ • s , ξ • s ) is a control equilibrium point of system (2), if there exists an input u s ∈ U such that The discrete control equilibrium sets X • s NL and X • s NL are the sets generated by every feasible control equilibrium is called a control equilibrium orbit of the NICS, and the set of all orbits VOLUME 10, 2022 associated to the discrete equilibriums sets (X One important property for control systems is the stability of the equilibrium points. In ICSs, the stability property is related to the orbit. Therefore, it is necessary to define this property as an extension of the well-know stability definition for dynamical systems: Following this definition, the next result establishes sufficient conditions for the asymptotic stability of the equilibrium orbit based on the asymptotic stability of the equilibriums of the underlying discrete systems. It reads as: Proposition 1: If a discrete control equilibrium ξ s of system (2) is asymptotically stable (under the standard definition for discrete systems), then the equilibrium orbit o s related to the NICS is asymptotically stable.
The proof of this result is detailed in the Appendix and in [42], but it relies on the assumption that f (ξ ) is Lipschitz, and hence, for a finite sampling time T s , the autonomous response is bounded.
For the implementation of the control strategy, first, consider an approximation of the nonlinear system (1) obtained by using the linear terms of the Taylor series of the system's dynamics f (ξ ) around the current state value ξ op at each time step. Note that point ξ op is not necessarily an equilibrium of f (ξ ), and thus,ξ ≈ẋ = f (ξ op ) + A(x − ξ op ). Therefore, the linearized system reads as: where matrix A = df (ξ ) dt | ξ op is the Jacobian of f evaluated at the current operation point, and matrix E = f (ξ op ) − Aξ op is the remainder of the first-order Taylor expansion. The linearized system (5) can be discretized to describe its behavior at times τ k and τ + k , obtaining a linear discrete-time subsystem of the form: . Note that matrices A • and E • must be updated at each time τ k since their continuous counterparts A and E depend on the current state ξ (t). It is also to remark that x • (k) = x(τ k ) describes the ICS before the input is applied, while is the state when the control input has already been applied. For this reason, despite both subsystems are necessary to characterize the ICS, the first subsystem x • is the one used in the control strategy to generate the control action [29]. In addition, as defined for the nonlinear case, the state x • s represents a control equilibrium of system (6) when satisfying All equilibrium points that satisfy the previous expression are associated to an equilibrium control set X • s . Finally, to represent the objective region where it is desired to steer the state ξ , the target set X Tar ⊆ X is defined. To accomplish the control goal, it is necessary to define a target counterpart of the equilibrium sets. In this regard, X •Tar s NL represents the maximal equilibrium set of subsystem (2) such that the equilibrium orbit set associated with all points in X •Tar s NL satisfy O Tar s ⊆ X Tar . For the discrete counterpart of the linearized system x • , a maximal equilibrium set X •Tar s associated with the target set X Tar is also defined.

C. CONTROL STRATEGIES FOR NICS 1) NONLINEAR MODEL PREDICTIVE CONTROL
In this Section, the ZMPC formulation to steer the state of a NICS to a target equilibrium set X •Tar s NL is described. The strategy is an extension of the linear impulsive case developed in [8] to the nonlinear case. The main characteristics of this strategy are (i) the capability of using a whole set as a target, (ii) the introduction of new decision variables x a , u a that are forced to be equilibriums of the system, and (iii) the update at each time step of matrices A and E used in the prediction model around the current operation state ξ op .
Given the current state ξ • (k), the optimization problem P MPC to be solved at each time step k is: where the first two terms of the cost function aim to steer the state and input to the artificial variables over the prediction horizon, and the last term of the cost function minimizes the distance between the artificial variables and the equilibrium target. In the prediction model, matrices A • and E • must be updated at each time k with the current operation state ), but they remain constant during the prediction horizon. The terminal constraint x • (k + H p ) = x a forces the state to reach the equilibrium at the end of the horizon, and the last constraint forces the artificial variables to be equilibriums of the linearized system. The solution of this problem is the optimal sequence U * = {u * (k), . . . , u * (k + H p − 1)} from which the first element u * (k) corresponds to the control law k MPC derived from the application of the receding horizon control policy, and it is applied to the NICS (1). A description of all parameters and variables introduced for the control strategy can be found in Table 1.
The following assumptions are established for the control strategy [8], [43]: 1) The NICS (1) is accessible, then the linearized model given by the pair (A, B) is controllable at each time step.
3) Matrices Q and R are positive definite, and P is a positive constant. 4) The target equilibrium set X •Tar s NL is such that it satisfies X •Tar s NL = X • s NL ∩ X Tar = ∅. Theorem 1 (Asymptotic Stability): Suppose that Assumptions 1) -4) are satisfied. Then, for any feasible initial state, the optimization problem (7) is recursively feasible for all time steps k and the control equilibrium set X •Tar s NL is asymptotically stable for system (1).
Proof: It has to be proven that the state ξ of the NICS (1) is steered and maintained in the target equilibrium set X •Tar s NL by means of the control law derived from the P MPC problem that uses the linear approximation of the system updated at each time step k. First, consider the change of variables x • = x • − ξ op and u = u − u op (where, for the impulsive case, u op = 0). Thus, the discrete linear system (6) used as prediction model can be rewritten as: In the definition of E • = ( Therefore, the discrete linear system results to be From this expression, the discrete system can be seen as a linear system with a deviation ( T s 0 e As ds)f (ξ op ), where the integral is a constant and f (ξ op ) is also known and bounded according to Assumption 2). Next, it is to consider that, when Assumptions 1) -4) hold, the optimization problem (7) has already been proven to be recursively feasible for linear and nonlinear systems in [8] and [43], respectively. Also, the set X Tar s (X •Tar s ) has also been proven to be an attractive generalized equilibrium set for the linear system (5) in [8]. Therefore, at each time step k, the ZMPC problem computes a control law K MPC that brings the state x closer to the target (due to the terminal cost and terminal constraint in P MPC ); hence, dist X •Tar (x(k)) → 0. By Assumption 3) and 4), dist X •Tar (x(k)) → 0 leads to f (ξ op (k)) → 0, then dist X •Tar s NL (ξ • (k)) → 0, and thus dist X Tar (ξ • (k)) → 0. Finally, since the discrete controlled system is asymptotically stable on X •Tar s NL , then, the NICS is asymptotically stable on X •Tar s NL , due to the asymptotic stability of the set of orbits O s associated to X •Tar s NL as stated in Proposition 1.

2) OFFSET-FREE NONLINEAR CONTROL FOR IMPULSIVE SYSTEMS
To improve the performance of the scheme in presence of parameter and modeling uncertainty, an offset-free feature is added to the ZMPC formulation. The offset-free strategy is based on the idea of correcting the prediction model and the target with a disturbance model that represents the plant-model mismatch such that z(k) → z t , when k → ∞, where z t is associated to a point ξ t ∈ X Tar . To this end, the system (1) is augmented with a model of a constant disturbance over time in the form: where d is a constant disturbance, and its corresponding matrices B d and C d represent the effect of the disturbance in the state and the output, respectively. These matrices are VOLUME 10, 2022 selected to guarantee the observability of the augmented impulsive system according to Corollary 2 in [29] and the observability of nonlinear impulsive systems in [44]. Afterwards, the linearized discrete-time subsystem in (6) and output equation are rewritten as: withx • = x • d • representing the augmented state, and the augmented matrices are given byÃ • = 0 e As dsB d . Since the augmented model (9) is also observable, the aim is to estimate both the state and the disturbance to correct the model in the control law. To that end, a state estimator (regardless of the strategy) is used and designed such that it is stable. Based on the characteristic of the systems here considered, in which there is a continuous nonlinear plant with measurements obtained at discrete points, it is decided, especially for applications, to use the hybrid Extended Kalman Filter (hEKF) [45]. The complete equations of the estimator can be seen in Appendix A, where the augmented model (8) is used to compute the a priori estimated state, and the augmented matricesÃ • ,C • are used to calculate the covariance of the estimated error and the filter gain. This estimator has been previously shown to converge to the state of nonlinear systems [45], [46].
It is to remark that, despite the constant dynamics of disturbance d, its value is updated each time k through the correction factor of the estimator, where the measurement information is introduced.
Next, given the current estimate of the augmented state, the optimization problem that solves the impulsive offset-free ZMPC (ZMPC-OF) scheme every time instant k is: where the dynamic constraint includes the correction of the prediction model using the augmented model, and the artificial reference is also corrected by including the estimated mismatch in the equilibrium constraint. Note that, by using this ZMPC formulation, both the prediction and the target are corrected in a single step, while when using standard MPC formulations, as the one presented for discrete-time systems in [12], it is required to solve two different optimization problems, one for the computation of the target, and then the MPC to compute the control action.
Theorem 2 (Offset-Free Control): Suppose that the MPC problem (10) is feasible for all k ∈ N, and the closed-loop system reaches an equilibrium ξ ∞ , d ∞ , z ∞ . In addition, consider the output target z t associated to ξ t ∈ X •Tar snl . Then, z ∞ → z t when k → ∞.
Proof: As it is assumed that the closed-loop reaches a steady state, the stability of the observer implies that the estimated state also reaches a steady state (ξ ∞ ,d ∞ ), which is associated to the pair ( Also, consider the target linear variables x t and y t = Cx t associated to the nonlinear counterparts ξ t and z t , respectively. Then, following the same approach as in [29], by defining δx ∞ =x ∞ − x t ∈ X , δu ∞ = u ∞ − u t ∈ U, and the offset = y ∞ − y t , from the equilibrium conditions of the estimator and the target it follows that (A • − I n x )δx ∞ + B • δu ∞ = 0 and Cδx ∞ = . In addition, by considering the change of variables δx(j) = x(j) − x t , δu(j) = u(j) − u t , δx a = x a − x t , and δu a = u a − u t , the optimization problem (10) can be rewritten to have the same form of (7) (see the complete development in [29]). As the MPC problem in (7) computes a stabilizing control law k MPC then, the only solution to the closed-loop system at steady state (A • − I + B • k MPC )δx ∞ = 0 is δx ∞ = 0, and hence, = 0, i.e., y reaches the target. Finally, as the state x is obtained at each time step around ξ , = 0 also implies that in the equilibrium z ∞ reaches z t .

THE ONCOLYTIC VIRUS THERAPY
This first application is conducted in the context of the experimental studies in [30], which investigated the antitumor efficacy of the genetically engineered oncolytic adenovirus ADPEDGHER against EG7 breast cancer cells in nude mice during 60 days. The mathematical model that describes the dynamics of OVT in [30] was developed in [35] and [36] and has the form: where S(t) is the number of susceptible tumor cells (10 6 cells), I (t) is the number of tumor infected cells (×10 6 cells), and V (t) is the number of virus particles (VP) (×10 9 virus) at the tumor site relative to the number of cells. These state variables are non-negative by definition. Model parameters have physiological meaning whose descriptions can be found in Table 2. The model (11) does not include immune responses because it represents experiments conducted in athymic nude mice exhibiting negligible immune responses [30]. The parameter = 0.001 is used to avoid singularity due to S(t) + I (t) = 0 whilst solving the system of ordinary differential equations.
The input u(t) of the model represents the injection of VP on predefined days τ k given by u(t) = u(τ k )δ(t − τ k ), where u(τ k ) is the amplitude of the control action, and δ denotes the Dirac function. The output of the model is the total number of tumor cells T = S + I , and is estimated from the tumor volume calculated using the relation 0.523 LeW 2 , where the length (Le) and width (W ) are measured using a caliper [30]. The density of the tumor is assumed to be 10 6 cells per mm 3 [35], [47].
To identify the model parameters and initial conditions, the burst rate of infected tumor cells was set as α = 3.5 to be consistent with previous data [35], [48], so the model (11) is globally identifiable. Afterwards, the parameter set (r, L, β, d I , d V ) was identified by fitting the model (11) to measurements of the tumor volume from 2 nude mice (denoted as subject 1 and 2) collected in [30]. The bestfitting parameters and initial conditions are given in Table 2. The goodness of fit statistics for each subject are R-squared (Subject 1: 0.95, Subject 2: 0.97) and Pearson correlation coefficient ρ (Subject 1: 0.97, Subject 2: 0.98) with p-value <0.001.
The accessibility of this system was analyzed in [6], where it was shown that the model (11)  The control objective of the OVT is to reduce (or eliminate) the total number of tumor cells T in less than 60 days with the lowest possible viral dose since toxicity is a major concern. However, in practical situations, physicians consider a successful oncolytic virus therapy if tumor cells are below 50 since, after this point, the tumor can be safely extirpated. To achieve the control objective, the oncolytic virus therapy is conducted using the proposed nonlinear schemes of estimation and control. Fig. 1 shows the state space evolution of ξ = [S, I , V ] under the nonlinear ZMPC scheme. From the state evolution (blue line), the characteristic behavior of an ICS can be visualized, where on predefined days a forced-response is produced, caused by the application of the viral doses, and then the free-response follows, where the tumor cells tend to grow given the absence of the control input. In addition, the state constraint set and target set are visualized along with the equilibrium sets of the nonlinear system and the approximated linear system every time step, i.e., each red line corresponds to the equilibrium of the linear system obtained at each time k. It can be seen how the linear equilibriums approach to the nonlinear equilibrium as the system evolves. This occurs thanks to the application of a control action computed to minimize the distance to the equilibrium inside the target at each time step. Thus, a point inside the nonlinear equilibrium X •Tar s NL is reached at the end of the simulation. Next, the tumor regression is compared under the linear ZMPC scheme for impulsive systems developed in [8], the nonlinear ZMPC in (7), and the ZMPC-OF in (10). For the linear formulation, a linearization of model (11) around a single operation point is required to be used as prediction model. The equilibrium point representing the healthy state (S = I = V = 0) would be desired, but as it is not controllable, a point close to it has been selected. For visualization purposes, S = 1.5 is selected, and the model equations are solved to find the equilibrium. As observed in Fig. 2, the linear ZMPC does not achieve the control objective for Subject 2, as the tumor does not reach the target zone. In contrast, when the nonlinear strategies are used, viral doses are increased at the beginning of the therapy, and the tumor decreases below 50×10 6 cells in less than 10 days, see Fig. 2. For subject 3, it is observed that both linear and nonlinear ZMPC strategies achieve the control objective. However, it is observed that with the linear ZMPC, the total number of tumor cells enters the target area a few days after it does with the nonlinear schemes. Therefore, the nonlinear ZMPC and ZMPC-OF allow oncolytic virus therapy to improve the control performance.
Since variations in biological rates may occur within a subject or between subjects, a sensitivity analysis was performed to determine the parameters that more affect the tumor cells' evolution and the direction in which they affect it (if they are positive or negative correlated) [49]. It resulted that the sensitivity of T tends to be larger for changes in r, β, d V , and α for subject 1, whereas the sensitivity of T tends to be larger for changes in r, β, d I , and d V for subject 2. The parameters r, d v , and L are positively correlated with variations in the total number of tumor cells, and the parameters β, d I , and α are negatively correlated with variations in the total number of tumor cells. The results support the need of personalized oncolytic virus therapy since subjects are affected differently by parameter variations.
Next, the results of the sensitivity analysis are used to assess the ZMPC strategy (7) in presence of parameter uncertainty to test the robustness of oncolytic virus therapy with that scheme. For each subject, the 4 parameters whose variations have the greatest impact on tumor cells were selected (i.e., those with the highest magnitude of sensitivity function), and varied by 10% and 20% from their nominal  values in the direction in which parameter variations increase tumor growth to simulate detrimental discrepancy. For instance, variations of +10% were considered for parameters r, d v , and L, while variations of −10% were set for β, d I , and α.
Even though it can be seen efforts made by the ZMPC strategy to achieve the objective, this strategy fails to enforce enough tumor regression in presence of given parametric variations, see Fig. 3. This occurs because the mismatches in biological rates are not taken into account in the prediction model nor the computation of the required dose for each injection, and thus, the ZMPC fails to adjust viral doses appropriately. These results motivate the design of strategies that compensate parameter uncertainty.
To overcome variations in biological rates, the formulation in (10) is applied to achieve zero offset. To this end, the matrices C d and B d of the augmented model are chosen such that the augmented system is observable. Then, the estimator (17) with the augmented matrices is designed to estimate the state variables S, I , V , and the disturbance d, and this information is provided to the ZMPC-OF to correct the mismatch between the dynamics in vivo and in silico.
Subsequently, the ZMPC-OF (10) was tested in presence of parameter variations. When parameter variations increase, the amount of virus administered tends to increase to compensate for the effect of the variations, see Fig. 3. Unlike the ZMPC strategy, the offset-free strategy enforces tumor regression in both subjects in presence of the simultaneous variations in the four most influential parameters. For subject 1, it is observed that the ZMPC-OF manages to steer the tumor cells into the target zone around days 9 and 12 when there are variations of 10% and 20%, respectively. For subject 3, it is possible to steer the tumor cells into the target around the fifth day in both variation scenarios.
The offset-free strategy adjusts viral doses dynamically to overcome variations in biological rates, see Fig. 3. These results suggest that oncolytic virus therapy benefits from the offset-free strategy to exhibit some robustness to uncertainty in biological rates.
Additionally, both control schemes (ZMPC and ZMPC-OF) were evaluated in the presence of modeling uncertainty in tumor dynamics. This is first tested by replacing the Gompertz growth function (logarithmic term) in the model (11) by rS(t), to consider that the tumor grows exponentially without treatment, as in [30] and [50]. The therapy guided by the schemes achieves tumor regression and the control objectives for all subjects in presence of the given mismatch in tumor growth dynamics, see Fig. 4. Although the ZMPC and ZMPC-OF reduce the tumor quickly with this mismatch, the ZMPC-OF delivered viral doses lower than those of the ZMPC, see Fig. 4. Since the ZMPC does not take into account the mismatch, the ZMPC output unnecessary high doses. In summary, the offset-free strategy is also advantageous for this mismatch in tumor growth dynamics since the ZMPC-OF reduces the doses while maintaining tumor regression, see Fig. 4. These results suggest that the ZMPC and ZMPC-OF allow oncolytic virus therapy to be robust to this mismatch in tumor growth dynamics.
The control schemes are also evaluated in the presence of stromal cells. Since stromal cells are few in tumor-bearing In Fig. 5, an offset is observed for subject 1 when controlling with the ZMPC in case 1, and a failure to reduce the tumor can be seen for case 2. In contrast, the ZMPC-OF achieved the control objective in both cases. It is also visualized that the ZMPC is more robust to model uncertainty in Subject 2 since the total number of tumor cells are steered to the target in both cases. Nevertheless, there is a risk of presenting an offset for higher variations in the model, as for case 2, tumor cells stabilize near the target boundary. On the other hand, when using the ZMPC-OF strategy, tumor cells are adequately reduced within the target zone, and this reduction is achieved even faster. These results suggest that the offset-free strategy is suitable for oncolytic virus therapy to exhibit some robustness to the detrimental impacts of stromal cells.

THE TYPE 1 DIABETES TREATMENT
In this second application, both control schemes are used to compute the insulin doses required to regulate the blood glucose (BG) concentration in subjects with T1D into the normoglycemia range defined as 70 ≤ BG ≤ 180.
To that end, the nonlinear model developed in [41] is used as prediction model and in-silico platform to perform the controllers. The model was modified from the original version to include the meal disturbance dynamics [51]. Thus, the state-space representation is: where x 1 represents the BG concentration (mg/dl), x 2 is the interstitial insulin (min −1 ), x 3 stands for insulin concentration in blood plasma (mU /l), and x 4 is the glycemia due to meal disturbances (mg/dl/min). The control input of the system u(t) corresponds to the insulin administered to the subject (mU /l/min). This input can be delivered as multiple daily injections or short duration pulses using an insulin infusion pump. In both cases, the input is more realistically modeled as an impulse, given the comparison with the sampling time.
Thus, it is given by is the magnitude of the insulin dose at time τ k . The second input of the model, r(t), accounts for food intake (mg/dl/min). The description of the parameters and the values used in this work can be seen in Table 3 [51], [52].
Regarding the accessibility of the system, it is to clarify that, as state x 4 only depends on the disturbance r, it is not accessible (there is no way to influence x 4 by moving u). Hence, the accessibility of model (12) is evaluated without considering this state. Next, it is to verify that the output h(ξ ) is not an autonomous element, and thus, that at least three doses are required to control the system. This is done by calculating the relative degree of the system, applying the condition (ii) of Proposition 2 in [1]: Based on this, it can be concluded that the impulse relative degree of model (12) is d 0 (y) = 3 and dL 2 f h(ξ ), g(ξ ) = 0 when x 1 = 0. Afterwards, based on Theorem 2 in [1], the accessible space can be determined as: whose determinant is: Therefore, the subsystem (x 1 , x 2 , x 3 ) is accessible if −p 2 3 x 1 = 0, and thus, x 1 and p 3 have to be different to zero. To show that the same condition holds to guarantee controllability of the linear system at each time step, the controllability matrix is calculated: whose determinant is p 2 3 x 1 , and hence, the same conditions as the accessibility are obtained. In the context of this application, these values are not zero, and therefore, the accessibility of the NICS and the controllability of (A, B) is guaranteed at all times.
The performance of the control strategies is assessed in a simulation scenario of 36 hours, where 3 meals of 6, 10, and 8 mg/dl/min are provided to the subject at 7:00h, 12:00h, and 19:00h, respectively. After the first 24h, no more meals are introduced to see the response of the system at steady-state (fasting condition). The target zone has been established as 90 ≤ BG ≤ 110, to steer the system into the normoglycemia zone without increasing the risk of hypoglycemia (BG < 70mg/dl). The sampling time has been set as T s = 5min, considering the treatment with an insulin infusion pump [38].
First, the evolution of glycemia and the insulin delivered is compared under the linear and nonlinear ZMPC scheme. For the linear ZMPC, the prediction model is obtained by linearizing in a single equilibrium point inside the target zone, while for the nonlinear ZMPC, the prediction model is updated at each time step as in Eq. (7). As observed in Fig. 6, the nonlinear ZMPC achieves a better performance as it manages to reduce glycemia levels at meal times faster than the linear ZMPC, and thus, a lower time in hyperglycemia is obtained. This shows the advantages of updating the model, since the prediction model is closer in each step to the actual behavior of the plant, and thus, a more accurate insulin dose can be delivered.
In Fig. 6, it can be visualized how the nonlinear ZMPC achieves to regulate the glycemia into the target zone when parametric uncertainty is not considered. But, in addition to the inherent errors in model identification, there are factors such as physical activity, the dawn phenomenon, and others, that alter glycemia behavior and thus alter the insulin requirements [10]. Therefore, the next step consists  of assessing the control strategy with parametric uncertainty. Here, parametric uncertainty of 30% is considered in the parameters related to insulin sensitivity (p 2 , p 3 ) and insulin degradation rate (p 4 ) to generate a plant-model mismatch. Each parameter is varied in the direction in which it provokes hypoglycemia. Then, it can be seen how the ZMPC fails to steer the BG levels into the target (there is a steady-state offset), and also, there is an episode of hypoglycemia after the postprandial time. This behavior is obtained since the ZMPC has no information regarding the mismatch in the prediction model, and despite it reduces the insulin doses with respect to the nominal case (given the feedback), it still cannot compute adequate doses to avoid the offset error. In contrast, from the same Figure, it is visualized how the ZMPC-OF reduces the insulin doses delivered to the subject to correct the differences detected with respect to the plant, and therefore, it achieves to drive BG levels into the target. Additionally, it is observed that the system has a better response during transient periods with the ZMPC-OF. This is visualized during meal times, when the strategy reduces the delivery of unnecessary high insulin doses and thus hypoglycemia is avoided.

III. CONCLUSION
This study presents a first approach to control schemes for nonlinear impulsive systems with application to relevant biomedical applications. The proposed control schemes have been formulated using an underlying discrete-time linear subsystem with matrices that approximate the nonlinear behavior at each time step. The first developed scheme is a ZMPC with artificial reference variables, with which it is ensured that at each time step the controller computes a control action that brings the state to an equilibrium inside the target. Next, the scheme was improved by adding an offset-free feature to better handle plant-model mismatches. This was achieved by augmenting the system with a disturbance model, then the extended state was estimated and this information was provided to the ZMPC to correct the prediction model and the target artificial variables.
In addition, this study shows how control engineering techniques bring sensible benefits to resolve clinical challenges in biomedical applications, where the administration of the drug can be seen as an impulsive input to the plant. In particular for oncolytic virus therapy, MPC has been shown to calculate personalized viral injections to achieve tumor regression in a therapeutic target zone. As expected, the robust performance of the ZMPC-OF is better than the one of the ZMPC. When variations in biological parameters or processes improve or impair tumor regression, the ZMPC-OF decreases or increases the control action consistently. This robust performance highlights trade-offs between viral doses and toxicity, which limit antitumor efficacy and robustness of therapies. As future work in oncolytic viral therapy, experimental tests with animals are envisioned, such as those carried out in [30], where tumor cells should be induced in mice and then treated with the virus doses calculated by the ZMPC or ZMPC-OF.
The benefits of the ZMPC-OF have also been shown in the context of T1D treatment, where the insulin doses are computed by taking into account the estimated mismatch between the model and the plant with the aim of steering glycemia to the desired zone. In addition, despite the objective of the ZMPC-OF is to deal with the steadystate offset, the control strategy has shown benefits to correct transitory responses, resulting in the avoidance of hypoglycemic events when the physiological variations in the patient tend towards this state. In the context of this application, the next steps of research are directed to clinical trials in subjects with type 1 diabetes, given that there are more technological advances for the implementation of a closed-loop system, where a continuous blood glucose sensor is used, and the dose calculated by the control strategy is administered through a continuous insulin infusion pump.
For further work in the development of more robust strategies to handle a plant-model mismatch that is constantly changing, it is possible to couple additional strategies to improve performance in the transients. For instance, both control formulations (ZMPC and ZMPC-OF) could be coupled with different adaptation methods: (i) direct adaptive methods to change control parameters such as the penalty matrices Q and R, or (ii) indirect adaptive methods to update the prediction model used in the MPC. For both methods, the changes could be based on the detection of variations in measurable variables, or historical data to detect patterns, in order to improve the performance indices appropriate to each application. In addition, with the offset-free MPC, there is an additional advantage consisting of the estimation of the plant-model mismatch, which represents an extra auxiliary variable that could be used to generate other adaptation laws or be coupled with fuzzy methods. Some works in literature that developed methods that could be coupled with the nonlinear ZMPC and ZMPC-OF are an adaptive fuzzy control in [53] and a sliding mode safety layer in [10].

APPENDIX A HYBRID EXTENDED KALMAN FILTER
The hEKF is expressed by: ξ − (t) = f (ξ (t)),ξ (t 0 ) =ξ + (k − 1) P − e (t) = AP − e (t) + P − e (t)A + Q e , P e (t 0 ) =P + e (k − 1) K e (k) = P − e (k)C (CP − e (k)C + R e ) −1 ξ + (k) =ξ − (k) + K e (k)[z(k) − Cξ − (k)] P + e (k) = (I − K e C)P − e (k)(I − K (k)C) +K (k)R e K (k). (17) where P e is the covariance of the estimation error, K e is the Kalman filter gain, and Q e , R e are the covariances of the process error and measurement error, respectively. Note that, the estimated state and its covariance P e are propagated between samples using the nonlinear dynamics of the system, and then, the estimation is corrected using the measured information at times k.

APPENDIX B PROOF OF PROPOSITION 1
This appendix shows the demonstration carried out in [42], on the stability of the orbits of a NICS when the discrete equilibrium points are stable. First, let us note that, if ξ s is an asymptotically attractive equilibrium of the system (2), then trajectories with initial conditions in a neighborhood of ξ s get closer as time advances. Based on this, it can be seen that in a neighborhood of ξ s , φ is a contractor, therefore: φ(ξ 1 , t) − φ(ξ 2 , t) < C ξ 1 − ξ 2 ∀t ∈ [0, T s ] (18) where the smallest constant C, known as Lipschitz constant is denoted as C φ .
After that, it is necessary to prove attractivity. Let > 0, as X * is attractive for the discrete system (2), then there exist ∈ X such that for every ξ • 0 ∈ there exist K = K (ξ • 0 , ) ∈ N such that dist X • s (ξ • (k)) < C φ ∀k ≥ K . Following the same argument shown before, it is possible to get that dist o s (o ξ • (k) ) < ∀k ≥ K . Therefore o s is attractive for the impulsive system (2).