# Physiological Control Law for Rotary Blood Pumps with Full-State Feedback Method
## 2 Materials and Methods
### Pulsatility Controller Design  
The dynamic equation that models the LVAD can be written as [9]  
\[\begin{array}{l}x(k+1)=A(k)x(k)+B(k)u(k)+\eta(k)\\ y(k)=C(k)x(k)\end{array}, \tag{9}\]  
where \(A(k)=\left[\begin{array}{cc}a(k)&0\\ d(k)&-c(k)\end{array}\right]B(k)=\left[\begin{array}{c}b(k)\\ 0\end{array}\right]\) and \(C(k)=\left[\begin{array}{cc}0&1\end{array}\right]\), here, \(x(k)=\left[\begin{array}{cc}x_{1}(k)&x_{2}(k)\end{array}\right]^{T}\), where \(x_{1}(k)\) is the pulsatility index of the pump speed (\(\overline{PI}_{\omega}\)); \(u(k)\) is the pump control input, which is the pulse width modulation (\(\overline{PWM}\)); \(\delta(k)\) is the system noise; and \(y(k)\) is the system output.  
The full-state feedback control is introduced the literature [28] in combination with the pulsatile flow estimator model to design the control law, \(u(k)=-Kx(k)\). The poles of the eigenvalues \(A(k)\) for the closed-loop system are modified to determine \((SI-\left(A(k)-KB(k)\right))\), where S is the Laplace variable. Therefore, the characteristic equation of the closed-loop system in Equation (9) is written as  
\[\det(SI-\left(A(k)-KB(k)\right))=0, \tag{10}\]  
\[\det(\left[\begin{array}{cc}S&0\\ 0&S\end{array}\right]-\left(\left[\begin{array}{cc}a(k)&0\\ d(k)&-c(k)\end{array}\right]-\left[\begin{array}{cc}K_{1}&K_{2}\end{array} \right]\left[\begin{array}{cc}b(k)\\ 0\end{array}\right]\right))=0. \tag{11}\]  
For a fast-dynamic response, we introduce the linear quadratic regulator optimal control [26]. Therefore, the control input is calculated in such a way that  
\[\rho=\frac{1}{2}\sum_{k=0}^{N}\left(x^{T}(k)Q_{1}x(k)+u^{T}(k)Q_{2}u(k)\right) \tag{12}\]is minimized subject to the control law \(u(k)=-Kx(k)\), where the \(Q_{1}\) and \(Q_{2}\) values are selected by the designer, in which their weighs are relative to the significance of the states and controls.  
Considering [PERSON]'s multipliers method in Equation (6), we can write the equation of the dynamic model of Equation (4) in terms of Equation (6), as follows  
\[\rho=\frac{1}{2}{\sum_{k=0}^{N}}\Big{(}x^{T}(k)Q_{1}x(k)+u^{T}(k)Q_{2}u(k)+ \gamma^{T}(-x(k+1)+A(k)x(k)+B(k)u(k))\Big{)} \tag{13}\]  
The partial derivative of this equation for three quantities of interest leads to  
\[\frac{\partial\rho}{\partial u}=u^{T}(k)Q_{2}+\gamma^{T}B(k), \tag{14}\]  
\[\frac{\partial\rho}{\partial u}=-x(k+1)+B(k)x(k)+B(k)u(k), \tag{15}\]  
\[\frac{\partial\rho}{\partial u}=x^{T}(k)Q_{1}-\gamma^{T}+\gamma A(k). \tag{16}\]  
Equations (14)-(16) represent the control, state, and adjoint, respectively, and the solution to these equations under certain boundary conditions results in the following:  
\[K=\Big{[}Q_{2}+B^{T}(k)S(k)B(k)\Big{]}^{-1}B^{T}(k)S(k)A(k) \tag{17}\]  
where S is calculated by solving the following Riccati equation:  
\[SA(k)+A^{T}(k)S-SB(k)Q_{2}{}^{-1}S+Q_{1}=0 \tag{18}\]  
This equation gives a quick convergence to the steady-state value; therefore, Equation (10) for the \(K\) value, representing the steady-state solution, is as follows:  
\[K_{\alpha}=\Big{[}Q_{2}+B^{T}(k)S_{\alpha}(k)B(k)\Big{]}^{-1}B^{T}(k)S_{\alpha }(k)A(k) \tag{19}\]  
where \(S_{\alpha}\) is the answer to the Riccati equation of discrete-time, and is given by  
\[S_{\alpha}=A^{T}(k)\Big{[}S_{\alpha}-S_{\alpha}B(k)Q_{2}B^{T}(k)S_{\alpha} \Big{]}A(k)+Q_{2}. \tag{20}\]  
This algorithm has been implemented as shown in Figure 3.