The resultant force is obtained by coupling the input and the servo forces. Both forces act on the rubber-material reaction disk, the deformation of which influences the component forces’ interaction (shown in Figure 2).

The relationship between the input and servo forces can be written by the following equation A 1 F servo + F Rd 2 + F k 2 − A 2 F i n − F Rd 1 − F k 1 = 0 (1) where A ₁, A ₂ are the cross-sectional area of the reaction disk’s center circle and outer ring. F _servo is the servo force provided by the motor. F _in is the force applied to the input rod by the driver. F _Rd1, F _Rd2 are the internal forces generated by the reaction disk’s deformation. F _k1, F _k2 are the spring force with the spring stiffness of k _pre . s ₁, s ₂, s are the displacement of the center-circle surface, the concentric surface, and the other side of the reaction disk as shown in Figure 2. Δ is the subtraction result between s ₁ and s ₂ (Δ = s ₁−s ₂).

The dynamic equation of TMC’s first piston can be written as x ¨ p = 1 m F servo + F i n − k m c x p − C x ̇ p − p m A m c (2) where x _p is the displacement of piston and is equal to the above-mentioned s. m is the piston’s mass; k _mc is the equivalent stiffness including the return spring and the tandem springs in TMC chambers. C is the coefficient of friction force. p _m is the master cylinder pressure. A _mc is the TMC cross-sectional area.

With the inlet valve as the border, a simplified hydraulic circuit can be divided into two chambers (shown in Figure 1). The first chamber indicates the closed volume from the TMC to the inlet valve, while the second chamber represents the wheel cylinder. The principle of equal-flow exchange exists between the chambers. The inlet and outlet valves are the actuators for regulating the flow. The inlet valves model can be expressed by q i v = C q λ i v A i v 1 ρ p m − p w i f p m − p w > − p crack − C q λ i v A i v + A c v 1 ρ p m − p w i f p m − p w ≤ − p crack (3) where q _iv is the flow rate passing the inlet valve and the direction from TMC to wheel cylinder is positive. C _q is the flow discharge coefficient. λ _iv is the orifice opening with a range from 0 to 1. A _iv , A _cv is the orifice passage area of the inlet and check valve. ρ is the brake fluid density. p _w is the wheel cylinder pressure. p _crack is the cracking pressure of the check valve. And the outlet valve is modelled by q o v = C q λ o v A o v 1 ρ p w (4) where q _ov is the flow rate passing the outlet valve; λ _ov is the orifice opening of the outlet valve with the value from 0 or 1. A _ov is the orifice area of the outlet valve.

The hydraulic model of TMC can be constructed based on the inlet valve’s flow rate. Its mathematical expression is as follows: p ̇ m = E m V m − A m c x p A m c x ̇ p − q i v p ̇ w = E w V w q i v − q o v (5) where V _m , V _w are the dead volume of the first and second chambers. E _m is the brake fluid bulk modulus. E _w is the equivalent bulk modulus of the wheel cylinder.

To facilitate the controller design, the mathematical model mentioned in this section discards some non-linear characteristics, including the following: i) the non-linear curve is replaced by a proportional straight line in the quasi-static increasing process; ii) there is a dead zone before the pressure is built, i.e., the stroke of master cylinder increases, but the pressure remains 0; iii) there is a hysteresis characteristic of the hydraulic system, i.e., the pressure of the piston in the same position on the forward stroke is greater than the pressure on the return stroke. As for the first item, as long as the linearization of the equiproportional ratio is applied appropriately, the effect on the actual control is relatively small. Meanwhile, the design of the pressure dead zone is mentioned later. Concerning the hydraulic hysteresis characteristics, the system stiffness in the return stroke is less than the forward stroke, which may cause the system to converge slowly in the return process.

H-shaped vehicles are more convenient for single-axle braking force control than X-shaped vehicles. The front-drive and rear-drive principles are the same for vehicles with single-axle independent control. In this paper, the front-wheel-drive (FWD) is considered as a research case, and its results can be extrapolated to single-axle-drive vehicles.

Hence, define the state variables x = x p x ̇ p F i n p m p w f p w r T and the input vector u = F servo q iv 1 q ov 1 q iv 2 q ov 2 T , the entire state space can be given by x ̇ 1 = x 2 x ̇ 2 = 1 m − k m c x 1 − C x 2 + x 3 − A m c x 4 + u 1 x ̇ 3 = g I x 1 , x 2 , u 1 x ̇ 4 = E m V m − A m c x 1 A m c x 2 − q iv 1 − q iv 2 x ̇ 5 = E w V w q iv 1 − q ov 1 x ̇ 6 = E w V w q iv 2 − q ov 2 (6) where g _I represents the mapping from x ₁, x ₂ and u ₁ to x ₃, which is determined by the state of the reaction disk. p _wf and p _wr are the front and rear wheel cylinder pressure. q _in1 and q _ov1 are front inlet and outlet valve flow. q _in2 and q _ov2 are front inlet and outlet valve flow. the signs of flow variables are positive when the flow directions come from the TMC to wheel cylinders.

The pressure in the TMC and wheel cylinder are decoupled during regenerative braking, thus the pedal feel is not always consistent with normal braking. The relationship between piston displacement x _p and the driver’s input force F _in is used to define the pedal feel in this study.

A MPC is developed to realize the precise regulation of both x _p and p _w , which indicates the decoupling of pedal feel and the wheel cylinder pressure. Another goal of the controller is to maintain brake-feel consistency during regenerative braking.

The non-linear model needs to be transformed into a local linear time-invariant (LTI) model by linearization and local simplification as the foundation of MPC design. Under ideal circumstances, the input force F _in is proportional to the piston displacement x _p . To ensure consistency of the pedal feel x ₃ can be written as follows: x ̇ 3 ≈ σ x ̇ 2 (7) where σ is the ideal pedal-feel coefficient between the displacement and input force.

When 7) is introduced into 6), by calculating the Jacobian matrix at the nominal point, the non-linear system can be converted to a linear system as follows: x ̇ t − f t 0 , x 0 , u 0 = ∂ f ∂ x t 0 , x 0 , u 0 x t − x 0 + ∂ f ∂ x t 0 , x 0 , u 0 u t − u 0 (8)

As shown in Figure 3, the hydraulic system is divided into three work modes during the regenerative braking process: fully hydraulic braking (FHB), distributed electro-hydraulic braking (DEHB), and coordinated electro-hydraulic braking (CEHB). The hydraulic system operates normally in FHB mode, with solenoid valves deactivated. The vehicle’s braking forces are all provided by the clamping force of the braking disc. When the motor’s regenerative braking is sufficient to meet the braking requirements of the driving wheels, the operating mode is switched to DEHB. The flow of the front wheels q _iv1 is limited to zero, and the rear outlet valve is used to regulate the rear wheel cylinder pressure. Similarly, in the CEHB mode, the continuously adjustable q _iv1 determines the front wheel cylinder pressure, allowing for blending braking of the driving wheels. It should be noted that the volume of a low-pressure accumulator is large enough that its internal pressure is almost zero.

The whole LTV plant is classified into three LTI sub-models in terms of control patterns. The state transition matrix set is solved by Eq. 8 as A k , B k , C k ⊆ A j , B j , C j j ∈ 1,2,3 , corresponding to FHB, DEHB, and CEHB in turn. The state at time k is manipulated by a state machine as shown in Figure 4.

Figure 5 illustrates a comparison of the linearized and original models. The relationship between x _p and F _servo has changed while q _iv is set to 0. The amplitude-frequency curve is moved down, and the phase frequency is also changed. Both linearized sub-models are capable of fitting actual curves.

The MPC controller is introduced to solve the non-linear problem. It’s a method for online optimization with a finite horizon. This is a multi-objective, multi-constraint optimization problem. MPC contains three features: model-based prediction, moving horizon optimization, and a framework that includes both feedforward and feedback. The non-linear electro-hydraulic model is simplified into the piecewise linear model in the previous section.

The controller’s objective is to implement multiple target tracking. The placement of MPC in the controller structure is shown in Figure 6. MPC receives the tacking reference from the upper-layer controllers, and then it’s the lower-layer actuators (EBooster and ESC) that conduct the commands from MPC. The control instructions are the increment signal to restrain the static error. The optimal objectives not only include the output references but also the input increments and input references. The increments decide whether MPC is conservative or aggressive. The input references can avoid the windup effect caused by the saturation of input and output.

The optimal constrained MPC problem can be formulated as follows: min z ∑ k = 0 N − 1 W y y t + k | t − y ref t 2 2 + W Δ u Δ u k 2 2 + W u u t + k | t − u ref t 2 2 s . t . x t + 1 = A j x t + B j u t y t = C j x t , j ∈ 1,2,3 Δ u min ≤ Δ u ≤ Δ u max Δ u k = 0 , k = N u , … , N − 1 u min ≤ u k ≤ u max , k = 0 , … , N u − 1 y min ≤ y k ≤ y max , k = 1 , … , N c (9) where W _y , W _Δu , W _u are the diagonal tunning weight matrices for output signals and input increments. N, N _u , N _c are the prediction horizon, the input horizon, and the constraints horizon, respectively.

According to the regenerative operating condition and control logic, a finite state machine is proposed to execute the state jump between the modes. The solution of MPC is presented in the single state where the state space can be regarded as the LTI in this following section.

Define x ̄ = x T y ref ( t ) T u ( t − 1 ) T u ref ( t ) T T and z = Δ u 0 ⋯ Δ u N − 1 T , the input sequence can be expressed by the input increment u 0 ⋮ u N − 1 = I u ⋮ ⋱ I u ⋯ I u z + 0 1 × x 0 1 × y ref I 1 × u 0 1 × u ref ⋮ ⋮ ⋮ ⋮ 0 1 × x 0 1 × y ref I 1 × u 0 1 × u ref x ̄ 0 (10) where the superscript denotes the dimension.

The output can be predicted by the state-space model as follows y 1 y 2 ⋮ y N = C C A B C ⋮ ⋮ ⋱ C A N − 1 B C A N − 2 B ⋯ C u 0 u 1 ⋮ u N − 1 + C A C A ⋮ C A N x 0 (11)

Then the simplification equations are obtained y 1 − y ref ⋮ y N − y ref = Φ 1 z + Γ 1 x ̄ 0 , u 1 − u ref ⋮ u N − u ref = Φ 2 z + Γ 2 x ̄ 0 (12) where Φ₁, Φ₂, Γ₁, Γ₂ are the constant matrices related to Eq. 10 and Eq. (11).

The weight matrices can be written by the weight vector Q ̄ y = b l k d i a g W y T W y ⋯ W y T W y ⏟ N (13a) Q ̄ Δ u = b l k d i a g W Δ u T W Δ u ⋯ W Δ u T W Δ u ⏟ N (13b) Q ̄ u = b l k d i a g W u T W u ⋯ W u T W u ⏟ N (13c)

The cost function is presented by J z , x ̄ t = 1 2 z T H z + x ̄ t T F T z + 1 2 x ̄ t T Y ⋅ x ̄ 0 (14a) H = 2 Φ 1 T Q ̄ y Φ 1 + Q ̄ Δ u + Φ 2 T Q ̄ u Φ 2 (14b) F T = 2 Γ 1 T Q ̄ y Φ 1 + Γ 2 T Q ̄ u Φ 2 (14c) Y = 2 Γ 1 T Q ̄ y Γ 1 + Γ 2 T Q ̄ u Γ 2 (14d)

When the constraints in Eq. 9 are transformed into the form of inequality, the standard form of QP is given by min z 1 2 z T H z + x ̄ F T z s . t . G z ≤ W + S ⋅ x ̄ (15)

This paper adopts an offline method to solve the QP. When the multi-parametric QP is from the MPC problem, the following matrix inequality is always satisfied Bemporad et al. (2002). Y F T F H ≥ 0 , H > 0 (16)

The optimal z is continuous and PWA linear functions of x within the critical region Bemporad et al. (2002). Based on the non-negative least squares (NNLS) solver Bemporad (2015), the local control law can be described by the PWA function. Δ u * t = F ̄ j x ̄ t + G ̄ j j : H ̄ j x ̄ t ≤ K ̄ j (17) where F ̄ , G ̄ , H ̄ , K ̄ is the approximate multiparametric solutions.

Mode CEHB is used as an example to draw a 2-dimensional projection of the polyhedral partition in the TMC piston displacement and piston velocity. The target displacement is 25 mm, and the remaining parameters without projection are set to 0. As shown in Figure 7, the polyhedra are divided into 240 blocks, and the actual trajectory converges from the initial point (0, 0) to the point (25, 0), via blocks 156 and 162 to block 1. This example shows that it is feasible to solve the QP problem by dividing the polyhedra offline.

The proposed MPC needs all states as the feedback signals. If a linear Kalman filter is used, due to the non-linear characteristics of hydraulic systems, the process noise produces cumulative errors with the dynamic process. It will lead to significant deviations in the observed values and bring disturbances to the MPC. If the disturbance is too large, the controller may not converge. In order to eliminate cumulative errors, the state-space matrice needs to be updated in real time. As a result, an extended Kalman filter (EKF) is adopted to estimate states in the non-linear transition process. Since the linearization has been mentioned in the previous section, it’s not described in this section. The extended recursive equations of EKF is as follows X ̂ k + 1 | k = A k X k | k + B k u k (18a) P k + 1 | k = A k P k | k A k T + Q k (18b) K k + 1 = P k + 1 | k C k T C k P k + 1 | k C k T + R − 1 (18c) X ̂ k + 1 | k + 1 = X ̂ k + 1 | k + K k + 1 y k − C k X ̂ k + 1 | k (18d) P k + 1 = I − K k + 1 C k P k + 1 | k (18e) where A _k , B _k , C _k are the value of Jacobian matrix for the non-linear transition process at the time k.

The measurement noise covariance matrix R is time-invariant matrix. There is a large difference in the process noise of different so that the process noise covariance is time-varying. This paper proposes an on-line approach for the process noise, which includes both the priori tunning value and Sage-Husa maximum posterior estimation algorithm Myers and Tapley (1976).

Once the regenerative mode occurs change, the process noise Q ̂ t is given the initial value Q 0 ( i ) where i denotes the work mode at time t. During the same operating mode, Q ̂ t can be updated by d k = 1 − b 1 − b k + 1 (19a) q ̂ k = 1 − d k − 1 q ̂ k − 1 + d k − 1 x ̂ k | k − A k x ̂ k | k − 1 + B k u k − 1 (19b) Q ̂ k = 1 − d k − 1 Q ̂ k − 1 + d k − 1 K k ε k ε k T K k T + P k − A k P k − 1 A k (19c) where k is the period number from time to the current. b is the forgetting factor and b ∈ (.95, 1). K _k is the current Kalman gain. ɛ _k is the residual as follows ε k = y k − C k x ̂ k | k − 1 (20)

Rubber material can be described by a non-linear GMM Jrad et al. (2012), which is divided into a quasi-static process and a dynamic process. Without the disturbance of rubber’s resilience, the static stressing and strain of rubber is shown in Figure 8. It’s the increasing curve with a decreasing slope, which can be fitted by a continuous three-piece linear function. By the means of the non-linear Trust-Region-Reflective Least Squares Moré and Sorensen (1983), the fitting piecewise function is obtained with the R ² of .9700. It includes 5 parameters with 2 x-coordinates and 3 slopes. The detailed solution process is not discussed in this paper.

The rubber’s dynamic characteristic is modelled by the Generalized Maxwell Model (GMM). A second-order transfer function with a zero point can be used to fit the reaction disk according to experimental data. Combined with the 3-piece linear fitting function, the frequency domain is presented in Figure 9. If the frequency is less than 10 Hz, the amplitude and phase can be regarded constant. Without considering the delay of controller, the phase frequency characteristics of the reaction disk reflect the hysteresis of the motor output force relative to the driver input force. When above 10 Hz frequency, the rubber hysteresis will exceed 20^◦. It means that with the higher frequency of the driver input force, the motor output force will obtain the greater delay. The gain of amplitude-frequency curves in different pieces are proportional and the phase-frequency characteristics are identical.

According to the fitting second-order model, the deformation difference Δ (introduced in Eq. 1) can be calculated by ξ ̇ 1 ξ ̇ 2 = A R d A R d ξ 1 ξ 2 + B Rd 1 B Rd 2 F i n F servo Δ = C Rd 1 − C Rd 2 ξ 1 ξ 2 (21) where ξ is the state of the second-order model; A _Rd and B _Rd are time-invariant; C _Rd is the function of input force and C R d ∈ { C R d piece 1 , C R d piece 2 , C R d piece 3 } .

The deformation difference can be measured by the displacement sensor. Its target value is acquired by the incremental updating method as Δ ̇ = C Rd1 A − C Rd 2 A ξ 1 ξ 2 + C Rd 1 B Rd 1 − C Rd 2 B Rd 2 F i n F servo (22)

The proposed MPC method can be used to solve another engineering problem. Before the pressure building are two non-linear processes: idle and dead zone [7]. Idle denotes that the mechanical part is not connected to the hydraulic part, while they are connected in the dead zone.

The non-linear process can be regarded as a specific condition of mode DEHB. The idle and dead zone can be converted to the rear wheels’ virtual flow, and the reference outlet valve flow is equivalent to the virtual flow simultaneously. The controller for the idle and dead zone is given by q vir = A m c x ̇ p u ̂ 3 = q vir m p c : 10 s . t . j = 2 0 ≤ x p ≤ x Idel + x Dead (23) where q _vir is the virtual flow of rear wheels; u ̂ 3 is the reference flow of rear outlet valves; x _Idel , x _Dead are the idle travel and dead zone displacement.