Sliding mode surface derivative using numerical method

In the article Controller Design using Bifurcation Map for Aircraft Spin Recovery, the sliding surface derivative is obtained using a numerical method.

The steady space system is defined by the following equations:

$$\dot{\underline{x}}=A(\underline{x},t)+B(\underline{x},t)\underline{U}, \\ \underline{x}=[V,\alpha,\beta,p,q,r,\phi,\theta], \\ s=\left[\begin{array}{c} \dfrac{d}{dt}(\phi-\phi_d)+\lambda_1(\phi-\phi_d) \\ \dfrac{d}{dt}(\alpha-\alpha_d)+\lambda_2(\alpha-\alpha_d) \\ \dfrac{d}{dt}(\beta-\beta_d)+\lambda_3(\beta-\beta_d) \\ \end{array}\right]$$

where $\lambda_1$, $\lambda_2$ and $\lambda_3$ are positive real numbers.

The article states:

The control law from reaching law in Eq. (6) can be derived as: $$\dot{s}=\dfrac{\partial{s}}{\partial{x}}\dot{x}=\dfrac{\partial{s}}{\partial{x}}(A+Bu)=-Qsgn(s)+Kf(s).\tag{8}$$

By solving the above equation (8) for u (3 x 1), one gets expressions for u as:

$$u=-\left(\dfrac{\partial{s}}{\partial{x}}B\right)^{-1}\left[\dfrac{\partial{s}}{\partial{x}}A+Qsgn(S)+Kf(s)\right].\tag{9}$$

The matrix $\frac{\partial{s}}{\partial{x}}$ is computed numerically. The control inputs computed are subject to position and rate constraints on control surface deflections as listed in the Table 1.

Which numerical method is used to compute $\dfrac{\partial s}{\partial x}$ in $u$?And also how would be the exact form of matrix of $\dfrac{\partial s}{\partial x}$ (element of $\dfrac{\partial s}{\partial x}$ in columns and raws)?

It's hard to say exactly which method they would use with certainty, but it's pretty reasonable to assume that they'd specify the method if it was anything special.

They easiest way to do this would be with a finite difference scheme. In this you approximate a derivative of a function by computing it at two locations that are a small, but known distance apart in your state space. The derivative (slope) is then the 'rise' (difference between the two values of the function) divided by the 'run' (the distance in state space)

The matrix in question is composed of the derivatives of each term of $s$ with respect to each variable in $x$. Treating $s$ as a function of $x$ you can readily form approximations of the derivatives as follows:

$$\frac{\partial s}{\partial \alpha} \approx \frac{s(V,\alpha +\Delta\alpha,\beta,p,q,r,\phi,\theta) - s(V,\alpha -\Delta\alpha,\beta,p,q,r,\phi,\theta)}{2\Delta\alpha}$$

Where I've introduced $\Delta\alpha$ as a small perturbation to give to nearby but different values at which to evaluate $s$. This will give will approximate s at $x = (V,\alpha,\beta,p,q,r,\phi,\theta)$ (between the two locations) and the error of the approximation will be proportional to ${\Delta\alpha}^2$.

• Dear Dan thank you for your answer.Useful.Could you please explain more about the ds/dx? Imean here there are 8 states and also 3 output signal to surface actuators .ds/dx matrix should be 3*8 matrix.I'm I right?(how would be the exact form of ds/dx related to the states (Raws and columns)?And also in S matrix we have for example [ d/dt(alpha-alphad)+lambda*(alpha-alphad)].Should we perturb alpha in both term (d/dt(alpha) & lambda*(alpha)? Commented Jan 8, 2016 at 8:37
• Yes, the $ds/dx$ matrix should be 3 rows, 8 columns. You can verify that from the multiplication by $B$ in the $u$ equation. Also, yes you should perturb $\alpha$ in both. Note you will not simply adding a 'bump' to the time derivative, rather you should compute the value at the new state space locations.
– Dan
Commented Jan 8, 2016 at 17:31
• Dear Dan thank you for helps.Something that is unclear to me is ds/dx elements.If we show ds/dx element s as [1 2 3 ...8;11 22 33...88;111 222 333 ...888] ، what whould be for example 3,33,333 according to your answer?Is it possible to add ds/dx formula (for each elements)to your answer? 1?11?111? ...888? Commented Jan 9, 2016 at 16:11
• Dear Dan as you know (ds/dx)*B be singular its inverse is none.Because of this I insist on the structure of ds/dx?thanks alot. Commented Jan 9, 2016 at 16:19