Fractional-Order Control of a Nonlinear Time-Delay System: Case Study in Oxygen Regulation in the Heart-Lung Machine

A fractional-order controller will be proposed to regulate the inlet oxygen into the heart-lung machine. An analytical approach will be explained to satisfy some requirements together with practical implementation of some restrictions for the first time. Primarily a nonlinear single-input single-output (SISO) time-delay model which was obtained previously in the literature is introduced for the oxygen generation process in the heart-lung machine system and we will complete it by adding some new states to control it. Thereafter, the system is linearized using the state feedback linearization approach to find a third-order time-delay dynamics. Consequently classical PID and fractional order PIλDμ controllers are gained to assess the quality of the proposed technique. A set of optimal parameters of those controllers are achieved through the genetic algorithm optimization procedure through minimizing a cost function. Our design method focuses on minimizing some famous performance criterions such as IAE, ISE, and ITSE. In the genetic algorithm, the controller parameters are chosen as a random population. The best relevant values are achieved by reducing the cost function. A time-domain simulation signifies the performance of PIλDμ controller with respect to a traditional optimized PID controller.


Introduction
During the complex heart surgery like heart transplantation it may be needed to temporarily interrupt the heart beat.As a result the heart gets cold and set to a resting condition.In this situation the blood, which provides enough oxygen, is not fed to the human body.However this quickly harms the tissues of the body.Even by applying a body-cooling technology (hypothermia), the tissue damage cannot be prevented for a bigger duration time.Thus an artificial machine is essential to take over the functions of the heart and the lung during the surgery.This machine is called heartlung machine (HLM) which will be located in extracorporeal circulation (ECC).The HLM receives the blood from vein, removes the carbon dioxide (CO 2 ) from the blood, and adds required oxygen to the blood whilst delivering to the artery.In 1960 the HLM machine was used in the heart surgery for the first time.It is now well improved and still is in use as a reliable and effective method [1,2].
An HLM is connected to the vascular system of the patient body by some specific tubes.During the cardiopulmonary bypass (CPB) surgeries the patient will be anaesthetized.Thus many vital and automatic tasks of the body will be interrupted.There is an alternative oxygenator as an artificial lung in the HLM machine to take over the function of gas exchange through propagation in a thin membrane with a high surface (1.8 m 2 ).
This gas exchange operations will be controlled by the gas fractions and the flow rate of the total mixed gases injected to the oxygenator.This operation is still regulated manually.Nowadays, there is no automatic commercial blood-gas machine.Since the gas exchange operation in the oxygenator is a complicated and nonlinear process with delays together with longtime constants, a manual control of this process is somehow difficult [1][2][3].This problem arises due to nonlinear and time-delay phenomena in the heartlung machine.Therefore an automatic approach must be taken in action for the gases regulation in the HLM.This is an auxiliary aim to investigate its performance.
PID controller is the most commonly used controller in industry due to its simple construction and regulation [4][5][6].This controller can be a useful tool in the heart-lung machine control.According to the development of fractional order calculus application in recent years, the researchers are trying to extend the common PID controller use to the fractional controller [7].This is almost because a complex model of the controller can be briefly introduced by a concise fractional model.A fractional order PID controller was first introduced by Podlubny [8,9] by defining it as PI λ D μ with λ, μ ∈ R.He has also shown that the stability and the performance characteristics of a system can be improved through choosing a real value for the integral and derivate orders.Frequency domain approach by using fractional order PID controllers was also studied in [10].Further research activities run in order to define new effective tuning techniques for nonintegerorder controllers by an extension of the classical control theory.To this respect, in [11] the extension of derivation and integration orders from integer to noninteger numbers provides a more flexible tuning strategy and therefore an easier achieving of control requirements with respect to classical controllers.In [12] an optimization method is proposed for tuning of FOPID controllers to achieve some specified desirable behavior of the controlled system such as gain margin and phase margin, high frequency noise rejection, and output disturbances attenuation.In this paper we use time-domain criterion such as integrated absolute error (IAE), integrated square error (ISE), and integrated time absolute error (ITAE).Consequently a fractional order controller is presented whilst the quality of the application in terms of the time response of heart-lung machine is also assessed.In this regard a genetic algorithm is used to find a set of the controller coefficient.This will be done through a fitness function to achieve rise in the performance indices of the rise time, overshoot, and settling time.
The rest of the paper is organized as follows.First a fractional order PI λ D μ controller is introduced.Then a regular model for the heart-lung machine is introduced involving time delays.
In the current literature several efforts have been made in the feedback linearization technique besides analytical approaches are provided to linearize and design appropriate controller hereafter.Therefore the model is reduced to a linear model by using a state feedback linearization technique.Then classical PID and fractional PI λ D μ controllers are applied to the linear system.Parameters of the PID are first tuned through Ziegler-Nichols method.Since the obtained parameters result a response with large overshoot together with long settling time, a better output response is of the request.Simultaneously, a best response of using different methods is compared in order to find a better outcome.Accordingly, three different criterions are used to achieve an optimum parameter set.Parameters of PID and PI λ D μ controllers are achieved using genetic algorithm in order to minimize the cost function.
Both PID and PI λ D μ controllers, where the gains are optimized by genetic algorithm, are applied on the heartlung machine.The quality of the system response is verified in a simulation study.

Fractional Order PI λ D μ Controller
A real value of the order of integral and derivative terms provides more flexibility in the controller design to achieve a desired response.This motivates many researchers to generalize the common controller to the fractional order one [9].PID controllers are the most widely used controllers in industry.A generalization of classical PID controller to the fractional order one is expressed as PI λ D μ by In fact an appropriate time domain action will be written as follows: The performance of PID and PI λ D μ controllers on the heartlung machine is investigated in the next section.

The Oxygenator Model
In this paper we will use the proposed model in [2].However the shortcoming of the previous model will be modified by adding the blood gas analyzer (BGA) dynamics.There is no time-delay seen in the achieved model.Thus it will be convenient for the linearization of the oxygen exchange process model.The linearization process results a system without delay time.Thereafter the produced delay by the mass transportation will be added to the open loop model as a pure delay.Ultimately PID and PI λ D μ controllers are designed for the time-delay linear system.
The proposed model in [2] explains the gas transportation process in the oxygenator.The model can be considered as a nonlinear single-input-single-output process where the blood flow rate q b may be treated as a predictable perturbation.However from the prior knowledge it is assumed constant at the beginning of the process.Here p CO 2 is neglected to be considered as a process output.In this paper, just the partial pressure of the blood oxygen is controlled.This system can be described as follows: where where u = pO 2,g,in = p bar FiO 2 is the control input.Meanwhile, the appeared vascular conditions apart from the blood flow rate q b , in the state variables are assumed constant.It is therefore assumed to be a perturbation that directly acts on the states.The state variables are considered as follows: The values of pCO 2 and pO 2 processes are considered as output.The equations in [4] are accordingly modified to: A full list of constant values is available in [2].The carbonic acid hydration is according to (18)-( 19) and the oxygen correction curve can be described through (20)-( 22): pO 2,virt = x x = x 6 10 0.024(37−Tb)+0.4(x3−7.4)+0.06log(40/x5) , ( 21) where a j , j = 1, . . ., 7, are constant parameters obtained from experimental observations.The output vector c in (3) delivers the measured by BGA oxygen partial pressure to the output.Without loss of generality the output is supposed as

Time Delay of the HLM System
The oxygen and neutral gasses mixing process in the mixer consists of a time delay T d1 (q g ) [2].This delay is caused by the metric distance between the mixer and oxygenator due to the gas flow rate in the tubing system.The blood gas analyzer (CDI500, Terumo, Japan) (Terumo Cardiovascular Systems Corporation, 6200 Jackson Road, Ann Arbor, MI 48103/800-262-3304/734-663-4145 (http://www.terumo-cvs.com/)).has an internal sampling time of T s,BGA = 6 sec with the sample blood 1 (mmHg).
The dominant time constant of T BGA is shown in (6).The next time delay is concerned with BGA which is indicated by T d2 (q b ) for to the blood flow rate.This delay is caused by the blood transportation from the artery, vein, and oxygenator to the BGA sensors.In this case no internal variables of the model are needed.The above two time delays can be integrated as T d = T d1 (q g ) + T d2 (q b ).

State Feedback linearization
The model in ( 6  model and by using the state feedback linearization method a linear model is achieved.Later, the time delay is added to the linear model.The linearizing control law using the state feedback is as follows [13,14]: where L f x 1 and L g L f x 1 denote Lie derivative which is as follows: (25) In fact the linearization procedure is to continue the output differentiation until the input term "u" appears in the d γ y/dt γ term such that the coefficient of "u", that is, is supposedto be nonzero.Thus, the resultant linear model is of order γ with "ν" as input and "y" as output: where γ is called the relative degree of the system.If the system in (26) has stable zero dynamics and γ as the relative degree of system then the state feedback law is It stabilizes the system exponentially with the following characteristic polynomial: The aim is to find the characteristic polynomial which is Hurwitz; that is, the coefficients are positive and the roots are all negative real valued.For the oxygenator system, three times differentiation is needed to make the input appears.
Thus the relative degree of the oxygenator system will be found γ = 3.For more details about the state feedback linearization method one may refer to [13,14].However the linear stable system is as follows: In which the parameters (β 0 , . . ., β 3 ) are chosen such that the input restriction and the sampling time T s,BGA = T d = 6 s is met.In other words it cannot be considered a linear system with a very fast response.The "best practice" rule that sample time should be 10 times per process time constant or faster (T s ≤ 0.1T p ) provides a powerful guideline for setting an upper limit on both control loop sample time and bump test data collection sample time [15].Therefore poles of linear systems are assumed to be located at s 1,2,3 = −0.015.
According to the feedback linearization technique a linear system with time delay is achieved.Consequently the linear transfer function is as follows: e −6s (s + 0.015) 3 . (30) In this system the total time delay isT d = 6 sec.The fraction input oxygen FiO 2 and the partial pressure of the oxygen pO 2 are considered as the input and output, respectively.

The Controller Design and the Simulation
A state feedback linearization of the oxygenator block diagram of the heart-lung machine is illustrated in Figure 1.
In this paper, classical PID and fractional PI λ D μ controllers are applied to the oxygenator system.One of the most commonly used methods in tuning of PID parameters is the Ziegler-Nichols method [16].In this method first the system is set in a closed loop using just with a pure proportional controller.The proportional gain is tuned such that the sinusoidal response is achieved at the output.This gain and the period of the sinusoidal response are treated as the critical gain (k cr ) and the period (T cr ), respectively.Then by using Table 1, PID parameters are, respectively, tuned.The appropriate output response using the Ziegler-Nichols method is stable; however, there are oscillations, overshoot, and rather long settling time at the output response (Figure 2).Thus, an optimal method should have been used in order to obtain a response with fewer settling time with lower overshoot.
Genetic algorithm is one of the optimization methods [17] which is based on the natural selection such as inheritance and mutation.In this algorithm variables which must be optimized are considered as gens.Accordingly set of all variables are defined as a chromosome or a person in the population.Each person or chromosome is considered as a point in the search space.This algorithm needs a cost function or fitness function in order to guide the procedure towards an optimum solution.During the classic PID controller tuning, each of parameters I, P, and D is assumed as gen.A set of genes establishes the chromosome or a person.Similarly for the fractional PI λ D μ controller adjustment parameters P, I, D, λ, and μ are individually assumed as genes and likewise for the chromosome.
In the genetic algorithm, first an initial population is randomly generated in the search space.Then some persons are randomly selected from the mating pool as parents to breed new children.The crossover operator is applied on their chromosomes during the breeding procedure.Fitness of each person is calculated according to the fitness function.A person/persons who has/have more fitness is/are considered for the next generation.The genetic algorithm uses the mutation operator in order to escape from local minimum.This operator provides the diversity of the population in the search space for each generation.After some iteration a near optimum solution will be achieved in the search area.In fact three different criterions, IAE, ISE, and ITAE are considered as optimization indexes or the fitness functions in order to find optimum parameters of PID and PI λ D μ controllers.These are defined as follows: where y d (t), y(t), and e(t) are the desired output (the reference input), the system output, and the error signal respectively.By minimizing these three performance indexes we seek which one of these criterions reduces not only the overshoot but also the steady-state error and the settling time.Finally the genetic algorithm yields optimum parameters which minimize the IAE, ISE, and ITAE indexes according to Table 1.

Simulation Results
Repeating the genetic algorithm with different initial population achieves an optimum parameters of PID and PI λ D μ controllers as in Table 1.

Investigation of the Time Indices of the System
The range of oxygen partial pressure is set to 100-200 mmHg during the surgery.In this paper the oxygen partial pressure is assumed 160 mmHg.Thus the reference input is considered as a step function with the amplitude of 160 mmHg.
Applying the classic and fractional PID controllers with optimum parameters in Table 1 generates the step response of the system as depicted in Figures 2, 3, 4, 5, and 6.To be more specific, the rise time, the settling time, the overshoot, and IAE values are also compared in Table 2. From Figure 2 it can be concluded that despite of the fast response in the Ziegler-Nichols technique there are oscillations and overshoot in the response and the settling time is also high.However, the obtained responses by the GA-tuned PID controllers with IAE, ISE, and ITAE criteria, are found rather better than the classic PID controller counterpart and the Ziegler-Nichols technique.Furthermore the overshoot is found smaller whilst the error is also converged to zero.As shown in Figure 6 the achieved response by the GA-tuned fractional PI λ D μ controller with IAE, ISE, and ITAE results in better response in comparison with the PID controller tuned by GA with IAE, ISE, and ITAE criterias in price of fewer rise and settling time.

Analysis of the Performance of the Controller in Presence of the Perturbation
The blood flow varies from 1 to 5 mL/sec during the surgery.Therefore it is considered as a step disturbance.For analyzing the performance of PID and PI λ D μ controllers, a step disturbance is applied to the system which is shown in Figure 7.The system response in presence of the disturbance is also depicted in Figure 8 whilst applying PID and PI λ D μ controllers.From Figure 8 it can be seen that the response which is obtained by GA-tuned PI λ D μ controller, converges faster to the steady state with lower overshoot.Additional improvements over the classic PID are also reported in Table 2.

Conclusion
In this paper a model of the heart-lung machine with the time delay is considered.The model is treated as a single-input single-output with time delay.The goal in this system is to control the oxygen partial pressure, using the ratio of the input oxygen FiO 2 .Primarily the system is linearized into a third-order system with time-delay.PID and PI λ D μ controllers are chosen to control the heart-lung machine.The Ziegler-Nichols technique is used to tune PID parameters.Applying this controller results an oscillatory response with large overshoot.Then the GA is used to optimally tune the PID and PI λ D μ controllers parameters to find a better response with lower overshoot and settling time by minimizing three different criterions IAE, ISE, and ITAE.The time indices are investigated when a step input and disturbance are applied distinctly.Simulation results are presented in Figures 2-8 and Tables 1 and 2. The results verify that PI λ D μ controller, tuned by minimizing IAE criteria in comparison with PID controller which minimized ITAE, provides better performance characteristics in terms of the transient and steady-state response.The outcome is also verified in presence of the step disturbance rejection.

Figure 2 :
Figure 2: Step responses of the oxygenator controlled by PID controller with Z-N tuning rule.

Figure 3 :Figure 4 :
Figure 3: Step responses of the oxygenator controlled by PID and PI λ D μ controllers with IAE criteria.

Figure 5 :
Figure 5: Step responses of the oxygenator controlled by PID and PI λ D μ controllers with ITAE criteria.

Figure 6 :
Figure 6: Step responses of the oxygenator with IAE, ISE, and ITAE criteria.

Figure 7 :
Figure 7: Blood flow as a step disturbance.

Figure 8 :
Figure 8: The system responses in presence of a step disturbance controlled by PID and PI λ D μ .

Table 1 :
Optimum setting of PID and PI λ D μ controller parameters.
p CO 2 , pH virt , [H] rbc , pCO 2,rbc , pO 2 , [HCO 3 ] rbc , . . ., pCO 2,pl , [H] pl , [HCO 3 ] pl , [HCO 3] indicates the bicarbonate concentration and [carb] is the carbamate concentration.pCO 2,g,out and pO 2,g,out stand for the partial pressure of CO 2 and O 2 in the gas phase, respectively."rbc" and "pl" subtitles denote the red blood cells and the blood plasma, respectively.The variables in (5) are replaced with the oxygenator equations.The first two equations state the BGA actions in terms of the state representation by ẋ1

Table 2 :
Indices of the system time response due to using PID and PI λ D μ controllers.