Optimal Force Allocation and Position Control of Hybrid Pneumatic–Electric Linear Actuators

Hybrid pneumatic–electric actuators (HPEAs) are redundant actuators that combine the large force, low bandwidth characteristics of pneumatic actuators with the large bandwidth, small force characteristics of electric actuators. It has been shown that HPEAs can provide both accurate position control and high inherent safety, due to their low mechanical impedance, making them a suitable choice for driving the joints of assistive, collaborative, and service robots. If these characteristics are mathematically modeled, input allocation techniques can improve the HPEA’s performance by distributing the required input (force or torque) between the redundant actuators in accordance with each actuator’s advantages and limitations. In this paper, after developing a model for a HPEA-driven system, three novel model-predictive control (MPC) approaches are designed that solve the position tracking and input allocation problem using convex optimization. MPC is utilized since the input allocation can be embedded within the motion controller design as a single optimization problem. A fourth approach based on conventional linear controllers is included as a comparison benchmark. The first MPC approach uses a model that includes the dynamics of the payload and pneumatics; and performs the motion control using a single loop. The latter methods simplify the MPC law by separating the position and pressure controllers. Although the linear controller was the most computationally efficient, it was inferior to the MPC-based controllers in position tracking and force allocation performance. The third MPC-based controller design demonstrated the best position tracking with RMSE of 46%, 20%, and 55% smaller than the other three approaches. It also demonstrated sufficient speed for real-time operation.


Introduction
With the recent growth in assistive robots, collaborative robots, service robots, and their applications, enhancing robot safety has become an important area of research. These robots need to satisfy higher safety standards than conventional industrial robots [1]. Various approaches have been investigated for raising their level of safety. These approaches include active methods that monitor the environment and robot states using machine vision and/or sensor fusion and try to avoid risky situations; and passive methods that try to develop inherently safe robots that are not able to cause high damage in cases of system failure or unpredicted situations. Using hybrid pneumatic-electric actuators (HPEA) to drive the robot's joints has been proposed as a solution to increase inherent safety without sacrificing precision [2,3].
Since the introduction of the HPEA idea in 1987 [4], different designs have been proposed to develop this type of actuator. They include a combination of a DC motor and a rotary pneumatic motor

Mathematical Model
In this section, mathematical equations are derived to model the system described in Section 2.1. The dynamics of the load mass are defined by: where L m , p x , p F , e F , and f F are the load mass, load position, the force from the pneumatic cylinder, the electric motor's force, and the total friction force, respectively. To keep the load motion range symmetric, its datum is located where the cylinder and the linear motor are in the middle of their range. The friction force is mainly from the pneumatic cylinder, and can be modeled as the sum of dry and viscous friction [19] Typically, the electric motor's current input to force output response time is considerably smaller than the controller sampling period and the response times of other components. Thus, the motor's internal dynamics can be neglected. If the relation between the motor force and the current input to the motor is precisely known, we can consider e F as an independent input to the position subsystem that can be provided accurately and instantaneously. The same assumptions cannot be made about the pneumatic force as it is a function of chamber pressures, which have significantly slower response time to the valve commands. In the next section, the pneumatic force's dynamic model is derived.
The pneumatic force is the result of the difference in the forces at the two sides of the piston, as indicated in (3). In this equation, P , A , and atm P are the absolute chamber pressure, effective piston area of the chamber, and the atmospheric pressure, respectively.
Deploying conservation of mass, conservation of energy, and the ideal gas law for each chamber leads to (4) and (5) [20]. In these equations, m  , T , and R are the mass flow into the chamber, absolute temperature of the air, and specific gas constant of air (287 J·kg −1 ·K −1 ), respectively. S L and

Mathematical Model
In this section, mathematical equations are derived to model the system described in Section 2.1. The dynamics of the load mass are defined by: where m L , x p , F p , F e , and F f are the load mass, load position, the force from the pneumatic cylinder, the electric motor's force, and the total friction force, respectively. To keep the load motion range symmetric, its datum is located where the cylinder and the linear motor are in the middle of their range. The friction force is mainly from the pneumatic cylinder, and can be modeled as the sum of dry and viscous friction [19] as follows: Typically, the electric motor's current input to force output response time is considerably smaller than the controller sampling period and the response times of other components. Thus, the motor's internal dynamics can be neglected. If the relation between the motor force and the current input to the motor is precisely known, we can consider F e as an independent input to the position subsystem that can be provided accurately and instantaneously. The same assumptions cannot be made about the pneumatic force as it is a function of chamber pressures, which have significantly slower response time to the valve commands. In the next section, the pneumatic force's dynamic model is derived.
The pneumatic force is the result of the difference in the forces at the two sides of the piston, as indicated in (3). In this equation, P, A, and P atm are the absolute chamber pressure, effective piston area of the chamber, and the atmospheric pressure, respectively.
Actuators 2020, 9, 86 4 of 15 Deploying conservation of mass, conservation of energy, and the ideal gas law for each chamber leads to (4) and (5) [20]. In these equations, . m, T, and R are the mass flow into the chamber, absolute temperature of the air, and specific gas constant of air (287 J·kg −1 ·K −1 ), respectively. L S and L DV are the stroke of the cylinder and equivalent length for the dead volume in each chamber, respectively. .
Using the model for compressible flow through the orifice area and the ideal gas assumptions, the mass flow rate through each valve can be modeled as shown in (6) [21].
Here, C f is the non-dimensional discharge coefficient of the valve, P u is the upstream pressure, P d is the downstream pressure, K is the ratio of specific heats for air (1.4), and A v is the valve orifice area. The valve orifice area changes as the spool moves. For Enfield LS-V05s pneumatic valves, it has been shown [20] that the relation between the valve orifice area and the normalized spool displacement (z) can be approximated by fitting a curve in the form of (7), where λ 1 , . . . , λ 7 are the coefficients calculated using experimental data points, and the least squares fitting method.
Since two-way valves are used in this study, each valve performs only one of the charging/discharging tasks; thus, 0 ≤ z ≤ 1. Regarding [20], this normalized spool displacement inside each valve can be modeled as: In (8), u v and τ v are the normalized command to the valve and the spool time constant, respectively. It should be noted that if a single three-way valve was used with each chamber, regarding (6) and (7), the pneumatic valve would switch between "charge" and "discharge" modes when the sign of spool displacement would change. Since the upstream and downstream pressures are different in each mode, there would be two distinct equations for the valve mass flow for charging and discharging. This would lead to a piecewise function for the mass flow rate of the valve. To be able to perform the controller in real-time, a linearized plant model with constant state space matrices is desirable (at least within the MPC prediction horizon). Therefore, a single mass flow equation is preferred to a piecewise function. Hence, the pneumatic circuit has been designed with a pair of two-way proportional valves for each chamber instead of a three-way proportional valve. With this structure (see Figure 1), the mass flow rates into the chambers are:

Controller Design
In this section, four different controllers are designed to address the position control and input allocation tasks for the HPEA-driven system. The first three controllers employ MPC for position control, while the fourth employs conventional linear controllers.

MPC1: Controller with Linearized Full Plant Model
Assuming that the full plant model of Section 2 and its parameters are known, an MPC law may be used for both position control and force allocation, outputting the electric force command and the inputs to the valves. The mass flow rates into the chambers are calculated using (6), which is a function of the upstream and downstream pressures. We term this approach MPC1. The controller structure is shown schematically in Figure 2.

Controller Design
In this section, four different controllers are designed to address the position control and input allocation tasks for the HPEA-driven system. The first three controllers employ MPC for position control, while the fourth employs conventional linear controllers.

MPC1: Controller with Linearized Full Plant Model
Assuming that the full plant model of Section 2 and its parameters are known, an MPC law may be used for both position control and force allocation, outputting the electric force command and the inputs to the valves. The mass flow rates into the chambers are calculated using (6), which is a function of the upstream and downstream pressures. We term this approach MPC1. The controller structure is shown schematically in Figure 2. The plant model must be linearized for use in the controller. The purpose of model linearization is to convert the nonlinear plant model into the linear discrete-time state space form of (10). The arguments k and k + 1 denote the kth and (k + 1)th samples.
To obtain the state space matrices (i.e., , , , Given (11), the continuous-time state space matrices can be calculated as shown in (12). The discrete-time matrices can be simply calculated by applying a discretization method to the matrices from (12). We used Euler's method to obtain the matrices as shown in (13). The plant model must be linearized for use in the controller. The purpose of model linearization is to convert the nonlinear plant model into the linear discrete-time state space form of (10). The arguments k and k + 1 denote the kth and (k + 1)th samples.
To obtain the state space matrices (i.e., A d , B d , C d , and D d ), Equations (1) and (4)-(6) are first used to derive the following vector equations: .
Given (11), the continuous-time state space matrices can be calculated as shown in (12). The discrete-time matrices can be simply calculated by applying a discretization method to the matrices from (12). We used Euler's method to obtain the matrices as shown in (13).
In (13a), I is the identity matrix and T s is the controller sampling period. The full matrices are presented in Appendix A. The subscript "op" indicates linearization about an operating point. This is a point that includes states and inputs, as shown in (14).
There are different options for selecting the operating point (14). One option is to use the current system states and input, i.e., and F e op = F e (k − 1). If this option is chosen, the operating point and the state space matrices need to be recalculated with every controller sampling period, passed to the MPC's optimization solver, and then, kept constant within the prediction horizon. Note that the mass flow rate through each valve is modeled using distinct equations for the choked and unchoked modes (see (6)). Having four valves and each valve having two modes means there are 16 distinct versions of f and g. With each iteration, the correct version of the functions must be chosen (depending on the valve modes) and used in (12) to update the state space matrices. This makes the resulting MPC too computationally complex to implement with common processors at the 10 ms or faster sampling period required for this application. To achieve a computationally feasible controller, we instead use a fixed operating point for the plant model linearization. This leads to only one version for each of the state space matrices (instead of the 16 versions with the dynamic operating point).
The following three assumptions are made to derive the linear model about a fixed operating point: (1) The states and inputs stay close enough to the operating point, (2) Valve spool time-constant (τ v ) is small enough compared to the controller sampling time that the spool's dynamics can be neglected and z = u v , and (3) Since the model is linearized about a constant operating point, the valve orifice area function in (7) may be replaced by a linear function of A v = λ 8 z + λ 9 so that this single point will convey more information about the trend of the whole function. The state space state, input, and output vectors are shown below.
From (3), it can be seen that the cylinder can have infinite combinations of chamber pressures that produce the same pneumatic force. The strategy behind defining the second output is to limit the pressure choices by keeping the average of chambers pressures near P mid = (P s + P atm )/2, which means requiring the second output to track a constant reference of r 2 = P mid − P Aop + P Bop /2. In this paper, the operating point elements have been picked as u op = 0, x p op = 0, and .
x p op = 0. The operating points for chamber pressures have been chosen such that they satisfy (P Aop − P atm )A A − (P Bop − P atm )A B = 0 and P Aop + P Bop /2 = P mid . Therefore, the pressure operating points are: The MPC1 law uses (10) to predict the future states. It solves the optimization problem (17) every sampling instant with a prediction horizon of N p points, by using a time-efficient convex optimization algorithm. In (17), k is the current sample number, r = x p,d r 2 T is the vector of the two reference inputs, and Q, R 1 , and R 2 are positive-definite gain matrices.

MPC2: Simplified Two-Loop Controller
The MPC1 approach must solve an optimization problem every sampling period, which is computationally expensive since it includes models of both the payload dynamics and pneumatic dynamics. As a faster alternative, with our second approach, termed MPC2, the MPC law neglects the dynamics of the pneumatic subsystem. In other words, the MPC2 law assumes the actual pneumatic force equals its desired value, i.e., F p = F p,d . This design requires a two-loop controller structure as the MPC2 block is responsible for controlling the position subsystem solely. An inner-loop pressure controller handles the pneumatic subsystem and defines the command signals for each valve. The schematic of this control system structure is illustrated in Figure 3.  MPC2 uses (1) to derive linear state space model as in (11) with the following states, inputs, and output. The corresponding state space model matrices are presented in Appendix A.
The optimization problem declaration is almost the same as MPC1, with the only difference being adding one more constraint (shown in (19)) to limit the changes in u(k) between two consecutive samples. This change is necessary as MPC2 ignores the inner dynamics of the pressure subsystems and its limitations.
The solution to these conditions is shown in (20). MPC2 uses (1) to derive linear state space model as in (11) with the following states, inputs, and output. The corresponding state space model matrices are presented in Appendix A.
The optimization problem declaration is almost the same as MPC1, with the only difference being adding one more constraint (shown in (19)) to limit the changes in u(k) between two consecutive samples. This change is necessary as MPC2 ignores the inner dynamics of the pressure subsystems and its limitations.
∆u min ≤ u(k + j) − u(k + j − 1) ≤ ∆u max , j = 1, . . . , N p Before feeding F p,d from the outer-loop position controller to the inner-loop pressure controller, the pressure distributor should define the desired chamber pressures. This is done by satisfying Actuators 2020, 9, 86 8 of 15 the two conditions (P A,d (k) − P atm )A A − (P B,d (k) − P atm )A B = F p,d (k) and P A,d (k) + P B,d (k) /2 = P mid . The solution to these conditions is shown in (20).
A conventional PI controller, with a bounded trapezoidal integral term, is used to control chamber pressures as shown in (21). Here, Γ denotes the integral bound.
where ϕ j (k) = K I,p The pressure controller in (21) returns a value for u j that needs to be normalized and distributed between the charging and discharging valves. To do so, the distributor function, Ω, is used as follows. where

MPC3: Modified Two-Loop Controller
MPC1 and MPC2 are two extremes in terms of dynamic models for the pneumatic subsystem. MPC1 uses a linearized model of the pneumatic dynamics, whereas MPC2 neglects these dynamics. The MPC law proposed in this section, termed MPC3, is a compromise between these two extremes that is intended to balance the trade-off between performance and controller computation time. The controller structure of Figure 3 with the independent pressure controller is adopted, but the MPC position controller is expanded to include a simplified model of the pneumatic dynamics.
To obtain the simplified model for predicting the pneumatic force, F p , we assume that the proportional term in (21) is dominant and that the derivative of each cylinder chamber pressure is approximately proportional to its associated valve command signal, i.e., .
Substituting (23) into (3) then gives the following simplified dynamic model for the pneumatic force: Finally, the discretized form of (24) is added to the state space model of MPC2 to produce the state space model used by MPC3. Its state vector is defined by (25). Its input vector and output are defined by (18a) and (18b), respectively. The corresponding state space model matrices are presented in Appendix A.

Linear Two-Loop Controller
As a comparison benchmark, a conventional linear position controller was also designed and simulated. The controller structure is the same as that of Figure 3, except the MPC2 block is replaced with a linear position controller. This controller structure is similar to the two-loop structure used in [3]. Two independent controllers are used with the pneumatic and the electric actuators, with the former using feedforward + PD controller and the latter using a PD controller. Since the electric actuator has a faster response but significantly lower force capacity, the feedforward terms that compensate for the modeled inertial and friction forces are fed to the pneumatic actuator controller, while the electric actuator's force capacity is saved for reducing the residual error. The equations for these controllers are shown in (26) and (27).

Results and Discussion
The four controllers were simulated using MATLAB m code and Yalmip [16] running on a Windows 10 PC with an Intel Core i5 processor. To challenge the tracking and force allocation performance, a relatively large payload mass of m L =10 kg was used. The pneumatic actuator parameters for an SMC CM2XB25-300 air cylinder [19] and an Enfield LS-V05s pneumatic valve [20] were employed. The simulated linear motor is an Aerotech BLMC-92. To prevent motor overheating, the electric actuator's force has been limited to 30 N, which is lower than the 44.5 N continuous force rating of this Aerotech motor. All controller parameters were manually tuned with the goal of minimizing the root-mean-square error (RMSE). This tuning employed a coarse search followed by a fine search to reduce the possibility of stopping at a local minima. Table 1 lists the parameters used. To keep the models within the MPC-based controllers in linear time-invariant form, the static and kinetic friction components are treated as unmodeled disturbances. Only the viscous friction is incorporated in the MPC models, for which a 20% uncertainty is included in the friction coefficientĈ v (i.e.,Ĉ v = 0.8C v ). Figures 4-7 show the results from each controller. They follow a reference trajectory made up of step inputs, a 1 Hz sine wave, then a faster 4 Hz sine wave. This trajectory was chosen to be challenging and led to multiple saturations of the valve commands and electric actuator force (at 1 and 30 N, respectively), as shown in Figures 4-7. For easier comparison, the tracking errors are plotted on the same axes in Figure 8.
Actuators 2020, 9, x FOR PEER REVIEW 11 of 16 To keep the models within the MPC-based controllers in linear time-invariant form, the static and kinetic friction components are treated as unmodeled disturbances. Only the viscous friction is incorporated in the MPC models, for which a 20% uncertainty is included in the friction coefficient ˆv ). Figures 4-7 show the results from each controller. They follow a reference trajectory made up of step inputs, a 1 Hz sine wave, then a faster 4 Hz sine wave. This trajectory was chosen to be challenging and led to multiple saturations of the valve commands and electric actuator force (at 1 and 30 N, respectively), as shown in Figures 4-7. For easier comparison, the tracking errors are plotted on the same axes in Figure 8.   To keep the models within the MPC-based controllers in linear time-invariant form, the static and kinetic friction components are treated as unmodeled disturbances. Only the viscous friction is incorporated in the MPC models, for which a 20% uncertainty is included in the friction coefficient ˆv C (i.e., ˆ0 .8 v ). Figures 4-7 show the results from each controller. They follow a reference trajectory made up of step inputs, a 1 Hz sine wave, then a faster 4 Hz sine wave. This trajectory was chosen to be challenging and led to multiple saturations of the valve commands and electric actuator force (at 1 and 30 N, respectively), as shown in Figures 4-7. For easier comparison, the tracking errors are plotted on the same axes in Figure 8.     From Figures 4-8, it can be seen that the linear controller has a longer settling time, bigger overshoot, and larger error peaks when subjected to the faster trajectories, such as 4 Hz. However, all MPC-based controllers undergo chattering in the actuator force at some points, while the linear controller produces the smoothest input force and valve commands, which may lead to longer valve and actuator lives in practice. Between the MPC-based controllers, MPC3 produced the most accurate position tracking. A quantitative comparison is presented in Table 2. Although the linear controller has the lowest accuracy of the three, it has a much faster calculation time. MPC1 requires the highest calculation time which is undesirable for real-time operation. MPC3 yields the best tracking performance with RMSE reduction of 46%, 20%, and 55% compared to MPC1, MPC2, and the linear controller, respectively. MPC3 also reduced the mean absolute electric actuator force by 9%, and the mean absolute pneumatic actuator force by 14% relative to the benchmark linear controller. These force reductions have the benefits of reducing the energy consumption and operating cost.   From Figures 4-8, it can be seen that the linear controller has a longer settling time, bigger overshoot, and larger error peaks when subjected to the faster trajectories, such as 4 Hz. However, all MPC-based controllers undergo chattering in the actuator force at some points, while the linear controller produces the smoothest input force and valve commands, which may lead to longer valve and actuator lives in practice. Between the MPC-based controllers, MPC3 produced the most accurate position tracking. A quantitative comparison is presented in Table 2. Although the linear controller has the lowest accuracy of the three, it has a much faster calculation time. MPC1 requires the highest calculation time which is undesirable for real-time operation. MPC3 yields the best tracking performance with RMSE reduction of 46%, 20%, and 55% compared to MPC1, MPC2, and the linear controller, respectively. MPC3 also reduced the mean absolute electric actuator force by 9%, and the mean absolute pneumatic actuator force by 14% relative to the benchmark linear controller. These force reductions have the benefits of reducing the energy consumption and operating cost. Actuators 2020, 9, x FOR PEER REVIEW 13 of 16

Conclusions
In this paper, three novel control algorithms for the position control of an HPEA-driven system are investigated. The first (MPC1) uses MPC with a linearized full plant model. The second (MPC2) uses an outer-loop MPC with an inner-loop PI pressure controller with a simplified linear model. The third (MPC3) uses the same structure as MPC2 with a modified model to find a balance in the tradeoff between the position controller performance and computation time. A fourth controller, employed as a benchmark, uses an outer-loop PD plus feedforward controller with an inner-loop PI pressure controller. MPC2 showed the best performance in reducing the input forces (and thus, the energy consumption and operating cost). However, modifying the model in MPC3 leads to the best position tracking performance. Although MPC2 and MPC3 are computationally intensive, they can still be run in real-time using commonly available off-the-shelf hardware.
This research will be continued by investigating the extension of the MPC-based position control and task allocation approaches to rotary HPEA actuators with nonlinear payload dynamics. To avoid the subjectivity of the manual tuning approach, we will also be investigating the application of global optimization strategies (such as particle swarm optimization) to obtain the controller parameters.

Conclusions
In this paper, three novel control algorithms for the position control of an HPEA-driven system are investigated. The first (MPC1) uses MPC with a linearized full plant model. The second (MPC2) uses an outer-loop MPC with an inner-loop PI pressure controller with a simplified linear model. The third (MPC3) uses the same structure as MPC2 with a modified model to find a balance in the trade-off between the position controller performance and computation time. A fourth controller, employed as a benchmark, uses an outer-loop PD plus feedforward controller with an inner-loop PI pressure controller. MPC2 showed the best performance in reducing the input forces (and thus, the energy consumption and operating cost). However, modifying the model in MPC3 leads to the best position tracking performance. Although MPC2 and MPC3 are computationally intensive, they can still be run in real-time using commonly available off-the-shelf hardware.
This research will be continued by investigating the extension of the MPC-based position control and task allocation approaches to rotary HPEA actuators with nonlinear payload dynamics. To avoid the subjectivity of the manual tuning approach, we will also be investigating the application of global optimization strategies (such as particle swarm optimization) to obtain the controller parameters.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Matrices for the State Space Models Used in MPC1, MPC2, and MPC3
The matrices for the discrete-time state space model used in MPC1 are as follows:  For MPC2, the discrete-time state space matrices are as follows: Finally, the discrete-time state space matrices used in MPC3 are as follows: