Design and validation of an MPC controller for CMG-based testbed

In the last years, Control Moment Gyros (CMGs) are widely used for high-speed attitude control, since they are able to generate larger torque compared to “classical” actuation systems, such as Reaction Wheels . This paper describes the attitude control problem of a spacecraft, using a Model Predictive Control method. The features of the considered linear MPC are: (i) a virtual reference, to guarantee input constraints satisfaction, and (ii) an integrator state as a servo compensator, to reduce the steady-state error. Moreover, the real-time implementability is investigated using an experimental testbed with four CMGs in pyramidal configuration, where the capability of attitude control and the optimization solver for embedded systems are focused on. The effectiveness and the performance of the control system are shown in both simulations and experiments.


Introduction
The growing interest in space exploration is focusing the research in the development of space architectures able to be independent and self-sustainable. Nowadays, space missions for satellites require the capability to rapidly change the attitude with high slew rate, both to ensure the success of the mission and to avoid catastrophic impacts among crafts (Yadlin 2011;Courie et al. 2018). In general the attitude control can be defined as the process aimed to stabilize and adjust the orientation of spacecraft during space maneuvers (National Academies of Sciences and Medicine 2016). Precise attitude tracking control of small satellites is a demanding task due to the limited hardware resources of sensors, actuators, and processors equipped on-board. During the years, attitude control systems based on Control Moment Gyros (CMGs) have been used for a large scale of satellites as they can generate a large torque compared to Reaction Wheels (RWs) systems. Though CMG systems can provide rapid slew capability and high pointing accuracy not using any of the limited propellant dedicated to the main propulsion system, an important difficulty is their inherent geometric singularity problem, regardless of how many CMGs the system is equipped with. Various efforts have been made to overcome this troublesome problem, and many CMG singularity avoidance methods have been developed, with different advantages and disadvantages (Kraft 1993;Wie 2005;Lee et al. 2007;Chakravorty 2012).
A configuration with CMGs and steering scheme for rapid maneuvering of Earth observation satellites is presented in paper (Yadlin 2011). The CMG singularity arises when the torque vectors of all the CMGs are aligned on a plane, so that the CMG system cannot produce a control torque in the normal direction of that plane. Different configurations of CMGs systems are compared in Miceli (2007), using a single PID controller with a steering law based on the Moore-Penrose inverse matrix to avoid the singularity condition. The analysis is particularly interesting for the comparison among different actuator system: the first two are CMG systems in pyramidal configuration with a different skew angles, the third system is a cluster of two CMGs and four RWs.
Even if Proportional Integral Derivative (PID)-based controllers are largely used in industrial applications, they are used in cases where the plant has a simple dynamics and where an appropriate proportional and damping gains are able to ensure a stable behavior of the control system. These controllers, combined with a Singular Direction Avoidance (SDA) steering law (Ford and Hall 2000), are usually designed for CMG-based space systems, as in Xin et al. (2015); Bedrossian et al. (2005), Ciavola et al. (2019). However, if the system to be controlled is more complex, the implementation of an optimal control strategy could be better. The optimal control solution is easy and simple to design, and the feedback gains can be computed by solving a Riccati equation. The most successful methodology among the optimal control environment is the Model Predictive Control (MPC).
As clear from literature (Wang and Boyd 2009;Mayne 2014;Mammarella et al. 2018), however, one of the main drawback of the MPC control scheme is related to the online optimization problem and to the difficulty of embedding a real-time solver for the on-board implementation. A practical solution, usually used for MPC implementation on embedded systems, is the offline evaluation of the control law, uploaded on-board by means of lookup tables Bemporad et al. (2002). In literature, gain scheduling design are usually proposed for linear time-varying (LTV) system, but for spacecraft control using CMGs there is a limitation due to the high number of independent time-varying variables. So, gain scheduling design can be unfeasible. Moreover, in space applications, to handle uncertainties and fast changes on the dynamics, this solution is usually not adopted. For this reason, the real-time implementation (translated in a real-time solver) is a key aspect for the on-board implementation of an MPC control system. MPC was first applied successfully for plants with slow dynamics (Qin and Badgwell 2003;Yu-Geng et al. 2013). However, in the last years, it was widely applied to different areas, as space applications. Indeed, only few control methodologies are able to address challenging tasks, required by space observations and missions, as deeply explained in Eren et al. (2017). The main advantages of using MPC controllers are related to two main aspects: (1) the MPC controller can directly calculate the angular velocity of the gimbals without calculating the external control torque, (2) input, state and hardware constraints can be included in the implementation. The first aspect enables us to avoid using the conventional steering laws with the singularity problem to drive the CMG system. Besides, the second aspect guarantees computational feasibility, secure operation, and the best possible performance of the system. As a consequence, the stability of the closed-loop system and its feasibility are guaranteed, while even the tracking performance is improved under the constraints.
In Guiggiani et al. (2015) a Model Predictive Controller for spacecraft attitude tracking with RWs actuators is proposed and the controller is designed for desaturation of the RWs, without fuel consumption, compensating the gravity gradient torques. In paper Ashok et al. (2016) the problem of attitude control of spacecraft using CMGs in pyramidal configuration is faced with the implementation of a method to avoid the singularity condition by means of Nonlinear MPC (NMPC). A singularity avoidance term is introduced in the cost function to be minimized, in order to generate a control input that moves away from the singularity when the CMG system approaches the neighbourhood of the singularity, and to rapidly converge to a desired attitude when the CMG system stays away from the singularity. Even if the implementation of an NMPC is treated, this paper still uses a steering law to calculate the value of the gimbal angular velocities. When a steering law is included in the closed-loop scheme, the controller generates the input and the steering law generates the gimbal rates value, avoiding the singularity condition, separately. So, it should happen that the generated torque differs from the controller torque, compromising the control performance in avoiding singularity. This problem can be resolved by generating directly from the controller the control input that does not achieve singularity.
In light of these circumstances, this paper proposes an MPC based attitude controller for the CMG system, which directly generates the gimbal rates so as to avoid the singularity problem possibly caused by the steering laws. The proposed method is based on an MPC method with additional features in Wada and Tsurushima (2016), which are virtual reference state and integrator state. Although the literature provides a linear MPC controller, it optimizes not only a sequence of control inputs, but also a virtual reference signal and an integrator state. These additional optimizations improve attitude tracking performance and reduce the steady-state error with the guarantee of recursive feasibility and constraints satisfaction. The further contribution of this paper is to verify the real-time implementability of the proposed MPC controller using an experimental testbed with four CMGs in pyramidal configuration. By floating the experimental equipment via spherical air bearing and by 1 3 matching the center of rotation and the center of gravity of the equipment, the testbed simulates spacecraft attitude dynamics. Since ground experiments of spacecraft attitude control are very difficult, and furthermore MPC controllers require implementation of real-time optimization, there have been a few references, which conduct experimental verification of MPC based spacecraft attitude control. We demonstrate both simulation and experimental results of the proposed method, and also show an optimization solver, which is easy to be translated and computationally efficient for on-board implementation.
The paper is organized as follows. In Sect. 2 the dynamics of the analyzed system is considered, focusing on attitude dynamics. The testbed system is described in Sect. 2.2. The MPC controller is introduced in Sect. 3, including the virtual reference and the integrator state. The real-time implementability of the proposed MPC controller is also discussed in Sect. 3.1. The Simulation results are described in Sect. 4 and the experimental test in Sect. 5. Finally, conclusions are drawn in Sect. 6.

Preliminaries
Some preliminaries on Linear Matrix Inequalities (LMIs) and on MPC control system are first introduced. The LMI approach is used for the definition of the controller gains and the Yalmip Matlab toolbox Lofberg (2004) is considered.
Given the symmetric matrices F i = F T i ∈ ℝ n×n , with i = 0, … , m , a LMI problem has the form where x ∈ ℝ m . The inequality symbol in Eq. (1) means that F(x) is positive definite, i.e. u T F(x)u > 0 for all nonzero u ∈ ℝ n . This LMI is equivalent to a set of n polynomial inequalities in x.
Moreover, it is possible to find nonstrict LMIs, with the form The LMI problem in Eq. (1) is a convex constraint on x and the set {x|F(x) > 0} is convex indeed. Multiple LMIs such as where denotes a diagonal matrix composed by the elements.
Whether the matrices F i are diagonal, the LMI F(x) > 0 is just a set of linear inequalities. Nonlinear convex inequalities are represented as an LMI problem using Schur complements.
where Q(x) = Q(x) T , R(x) = R(x) T and S(x) depend affinely on x, is equivalent to It is often possible to deal with problems in which the variables are matrices as follows where A ∈ ℝ n×n is given and P = P T is the variable. We can put the Lyapunov inequality in Eq. (5) in the form of Eq. (1) considering P 1 , … , P m as a basis for symmetric n × n matrices (m = n(n + 1)∕2) and then taking F 0 = 0 and F i = −A T P i − P i A . Furthermore, leaving LMI inequalities in a condensed form as Eq. (5) may lead to more efficient computation. Consider the quadratic matrix inequality where A, B, Q = Q T , R = R T > 0 are given matrices of appropriate sizes and P = P T is the variable. It can be expressed as the LMI problem This representation shows that the quadratic matrix inequality in Eq. (7) is convex in P. In some problems it is possible to find linear equality constraints on the variables such as where P ∈ ℝ k×k is the variable. In order to write Eq. (8) in the form F(x) > 0 we can eliminate the equality constraint considering P 1 , … , P m as a basis for symmetric k × k matrices with trace zero (m = (k(k + 1)∕2) − 1) and P 0 as a symmetric k × k matrix with TrP 0 = 1 . Thus,

Model definition
As briefly introduced before, the main objective of this paper is the precise attitude control of a testbed, equipped with four Control Moment Gyros (CMGs) in pyramidal configuration. The testbed is considered as a rigid body, in which Euler's equation and quaternion kinematics are included Markley and Crassidis (2014). After the definition of the dynamic and kinematic equations of the testbed, the Linear Time-Invariant (LTI) mathematical model is derived in a state space formulation.

Attitude dynamics modeling
For the study of spacecraft dynamics, two reference frames are usually considered: an Inertial Reference Frame and a Body Reference Frame, as in Fig. 1.
The Inertial Reference Frame considered in this study is the Earth-Centered Inertial (ECI) that is also called Inertial Geocentric Reference Frame (Fig. 1). It is commonly used to study the motion of a body orbiting around Earth and it has its origin in the Center-of-Mass (CoM) of Earth.
The Body Reference Frame (Fig. 2) has been adopted to describe the attitude dynamics once introduced the CMGs to the equations of motion. The origin of the Body Reference Frame is in the CoM of the spacecraft and it is rigidly attached to the body.
Taking into account the previous assumption, the Euler's equation of motion can be written as Srinivasan et al.  where the subscript B indicates variables in the Body Reference Frame, B ∈ ℝ 3 is the angular velocity of the spacecraft, h B ∈ ℝ 3 is the angular momentum, M B ∈ ℝ 3 is the torque vector defined as F i is the force which has a momentum with respect to the CoM, and r i is the distance between the CoM and the considered force. For a generic body, the inertia tensor J directly relates the angular velocity with the angular momentum. The angular momentum h of a rigid body, referred to its CoM, can be defined Considering the testbed system, Eq. (10) can be rewritten as Thus, the total angular momentum is the sum of two terms: (1) the angular momentum of the spacecraft, and (2) h C B is the term given by CMG system, which for a pyramidal configuration is with where i ∈ ℝ 4 is the gimbal angle, for i = 1, ..., 4 ., and h = J ∈ ℝ denotes the angular momentum of each wheel. Since we suppose four identical wheels rotating at the same constant speed, h is a constant scalar.
Combining the previous equations, the Euler's equation here considered is For a cluster of 4 CMGs, the internal angular momentum vector h C B is a non linear function of the gimbal angles i , with i = 1, ..., 4 , of each CMG, as in Eq. (12). The time derivative of Eq. (12) is where the matrix A( ) ∈ ℝ 3×4 , considering a fixed skew angle , is The columns of matrix A( ) express the contribution of each CMG to ̇h C B . Finally, the following assumptions are considered 1. The initial angular momentum is zero, and consequently B (0) = 0 2. The initial torque applied to the body is zero, M B (0) = 0 With these hypotheses, starting from Eq. (15) and substituting Eq. (16), we can easily obtain the following simplified equation for the angular acceleration This latter equation relates the angular acceleration of the testbed to the angular rate of the gimbals ̇ , which is the control input vector in this system. For the definition of the model, the subscript that indicates the Body Reference Frame is omitted for simplicity.
The kinematic equation of the system using a quaternion representation is is the vectorial part of the quaternion, and q 4 is the scalar term.
The matrix (q) ∈ ℝ 4×3 is defined as: Combining Eqs. (18) and (19) obtained previously, the following mathematical model can be obtained where the system state vector is defined x p = q ∈ ℝ n p , with n p = 7.
Rearranging the previous equation, the state equation of the system is The state of the system is x p ∈ ℝ 7 and the control input is ̇∈ ℝ 4 , that is the angular rates of the gimbals. The matrix B p can be considered to be a (external) parameterdependent matrix as B ( ).
A linearisation is performed assuming that the initial gimbal angles are where y = [q T v , T ] T ∈ ℝ 6 and I 3 ∈ ℝ 3×3 is the identity matrix.

Testbed description
The testbed considered in this work is the one developed in Higashiyama et al. (2020), which is depicted in Fig. 3. Physical parameters of the testbed is summarized in Table 1. The testbed simulates a spacecraft using a metal plate equipped with the four CMGs in pyramidal configuration, a PC, two microcomputers, and three sets of counter weights. The metal plate is floated by compressed air via spherical air bearing, and the center of rotation and the center of gravity of the testbed is matched by manually adjusting the counter weights. The attitude and the angular velocity are estimated via a triaxial gyroscope and accelerometer called the 3-Space Sensor made by YEI Technology, and a low-pass filter. Moreover, the nominal moment of inertia is estimated by applying the least-square method to a regressor representation of the Euler's equations of the uncontrolled rotational motion. (21)

Design of an MPC controller
Starting from the research of Wada and Tsurushima (2016), a linear MPC with integrator as a servo compensator is considered to control the attitude of a spacecraft/ testbed. The integrator is included to reduce the steady-state error. From the optimization procedure, instead of only the control input, a virtual reference signal and an integrator state are also evaluated. The virtual reference signal w is a variable that substitutes the original reference signal r in the control algorithm in order to guarantee input constraints satisfaction. The controller has to decide at each sample time if the integrator state x c must be optimized or not. If x c is not optimized, it evolves according to the integrator equations in order to guarantee the steady-state error elimination. Instead, if x c is included in the optimization problem, the optimized value of this variable gives a rapid decrease of the cost function value. The plant system is described in the following state space formulation  where x p ∈ ℝ n p is the state of the plant, u ∈ ℝ m is the control input and y ∈ ℝ n y is the system output. In our case the matrix D p is null. In the MPC control system, the matrices of Eq. 21 are translated in discrete time as follows A p ∶= I d +Ā p dt , and B p ∶=B p dt . C p = C p . Moreover, following the literature, (⋅)(m|k) means a value at time m calculated at time k.
In our MPC formulation, both input and state constraints are considered as follows: where the matrices u ∈ ℝ 2m×m , u ∈ ℝ 2m , x ∈ ℝ 2n p and x p ∈ ℝ 2n p ×n p .
Finally, the integrator equations are with e ∈ ℝ n y . Combining the equations of the plant and the integrator, the error dynamics is defined as The matrices and are calculated for the tracking control problem in case of linear systems, and to guarantee input and state constraints satisfaction. The calculation procedure and the proofs are presented in Wada and Tsurushima (2016).
In this paper, the following cost function is considered and implemented where H s denotes the prediction horizon, and for a generic vector x and a square matrix A of appropriate dimension, ||.|| A denotes the weighted Euclidean norm, namely, ||x|| A ∶= (x T Ax) 1∕2 . The matrices Q, R, M, and P in (34) are positivedefinite. The first three matrices are design matrices, while the matrix P has to be obtained solving an LMI problem for the terminal controller of the dual-mode MPC guaranteeing recursive feasibility and stability, as described in Wada and Tsurushima (2016). The novelty in the defined cost function is found out in the last term, || (w(k|k) − r(k|k))|| 2 M to include the virtual reference signal w into the optimization procedure. The virtual reference changes from the original reference r in order to satisfy the constraints and to obtain better control performance. The last term prevents w from being unnecessarily away from r, and also guarantees that the deviation between w and r remains finite. Under the regulator condition relating to the existence of the matrices and , the convergence of w to a step reference r is shown in Wada and Tsurushima (2016).
Using the prediction matrices S and T, the state prediction of the error system up to the horizon H s is defined as wherẽ is the vector of the future state, instead ̃ is the optimized vector of the error system inputs, obtained by the optimization procedure. Now it is necessary to rewrite the constraints in function of the variables of the error system. The input constraints must be rewritten, considering The state constraints are rewritten as As already said, in the considered MPC, the integrator state can be or not optimized. In order to decide if x c has to be optimized or not, a threshold s must be (34) defined. The integrator state x c is reset at each sampling time so that the cost function value J(k) is minimized while J(k) is greater then s . The control algorithm with integral action is used to achieve zero steady-state error when J(k) is less than or equal to s .
The value of s is determined by trial and error procedure. If s = 0 , the robustness against the disturbances cannot be assured. However, if s is too large, the robustness against disturbances is improved but the tracking performance are not good. So, the parameter s has to be chosen by considering a trade-off between tracking performance and robustness against disturbances.

Real-time implementability
For the implementation of the MPC control system, the required resources are correlated with the problem size (i.e., number of constraints, prediction horizon) and, at the same time, the efficiency of the solver algorithm for the optimization of the cost function. In addition to these two computational indices, on-board hardware capabilities and software architectures should also be considered. Implementation issues are discussed in this Section due to hardware constraints and limited computational power of the embedded controller. Some simplifications of the code and analysis of the solvers are required. Moreover, since the GNC hardware is tested on Matlab/Simulink environment, the selection criteria for the solver includes compatibility with this environment. In this paper, in implementing the proposed controller to the testbed, we equivalently convert the optimization problem formulated as an LMI optimization problem in Wada and Tsurushima (2016) into a quadratic programming problem, which has much less computational burden, as where (k) is the optimized vector.
Moreover, the cost function is re-written using Matlab notation, to simplify the implementations on this environment. So, we have where (k) is a vector composed by the measured state of the system and by the reference signals.
The tested solvers are • quadprog, the interior-point-convex algorithm provided by the Matlab Optimization Toolbox to solve quadratic programming problem, • fmincon, Matlab function that can implement four different algorithms: interiorpoint, Sequential Quadratic Programming (SQP), active set, and thrust-regionreflective. Also this solver deals quadratic programming problems.
• quadwright, a quadratic programming solver developed by J. Currie presented in Currie (2014); Wright (1996), able to speed up the computational capabilities for embedded applications.
All these solvers are tested in simulations. The software quadprog can use a thrustregion-reflective method for handling problems with bounds and linear equality constraints. However, it showed an high computational effort and cannot be translated in C code. The solver fmincon is usually used for nonlinear constrained problems and it finds a constrained minimum of a scalar function of several variables starting at an initial estimate. This solution is not a global minimum, but it is computational efficient. Finally, the quadwright provides the smallest computational times, because it has been developed with a focus on efficient memory use and high speed convergence. For the experimental verification, we create a C++ code with Eigen and eiquadprog libraries. The Eigen library is a C++ template library for linear algebra, and the eiquadprog Guennebaud et al. (2011) library is a C++ routine working with Eigen data structures for solving quadratic programming problems based on the Goldfarb-Idnani active-set dual method.

Simulation results
The simulations are performed in Matlab, with an Intel Core i7 at 2.4GHz and 8GB RAM. The selected sampling time is 25 ms and the prediction horizon H s is 5 steps. The MPC solver, as previously said, is quadwright. The controller parameters are defined in Table 2. The reference signals for the quaternion are evaluated from the evolution of the angular velocity, starting from Eq. (19), for all different simulation cases. Furthermore, in all the simulations the initial conditions are null states x, except for x 4 (0) = q 4 = 1 and null gimbal angles . So, the initial conditions are the ideal quaternion and null angular velocities. Two set of simulations are considered: (1) in the first set of simulations, all the three angular velocities are changed. Trapezoidal variations of the angular velocities are simulated, (2) in the second set of simulations, triangular variations of an angular velocity is examined.   Figure 8. The control input of the MPC with integrator state and virtual reference is shown in Figure 9. In a similar way, noisy variation of the gimbal angles can be observed in Figure 10, whereas no noisy response is observed in Figure 11. Instead, no difference can be observed in the quaternion behavior, as in Figures 4, 5 and 6. Finally, to show the effectiveness of the proposed controller, a singularity index is evaluated, which is det(AA T ) , as in Higashiyama et al. (2020). Highest is the singularity index, more effective is the singularity avoidance. The singularity indices are in Figures 12 and 13. As in the previous results, a noisy behavior is observed in Figure 12. Moreover, a small value is obtained if no additional optimizations are included.
In the second simulation, as before, a step variation of the three angular velocities is considered with two variations of the angular velocities. A simulation time of 10 seconds is considered: (i) d = [2, 2, 2] T deg/s between 2 and 4 seconds, (ii) d = [4, 4, 4] T deg/s between 5.5 and 8 seconds. As for the previous case, the state variables (angular velocities and quaternion) are analyzed for the classical MPC and the one with optimization of virtual reference signal and integrator state. Moreover, control input and gimbal angles are also considered. All the results  Figs. 14,15,16,17,18,19,20 21. Finally, the singularity index is analyzed to show the effectiveness of the proposed MPC controller. As in Figs. 22 and 23, even in this simulation the singularity index of the MPC with optimization of virtual reference signal and integrator state is greater than the classical MPC with no additional optimizations.
These simulations demonstrate the effectiveness of the MPC controller with additional features, that optimize the virtual reference signal and the integrator state.

Second set of simulations
In this set of simulations, triangular variation of one angular velocity is proposed. This set of simulations is performed in simulation and by experimental tests as well (as detailed in Sect. 5). As in the previous section we show that the MPC with integrator state and virtual reference has better performance than the classical one. Thus,  As from Figs. 24 and 25, the controller is able to track the desired behavior, with a small control authority (Fig. 26). Moreover, the singularity index (Fig. 27) shows the ability of the controller to avoid singularity.
Similar results are obtained for the rotation around Z-axis (Figs. 28,29,30 and 31). The singularity index shown in Fig. 31 implies that if the value of the index is close to zero, the CMG system is close to a singular state. Thus, the higher is this index, the higher is the controller ability to avoid singularity.

Experimental results
After performing simulations in Matlab environment, the MPC with optimization of the virtual reference signal is converted in C++ program language, and applied to the testbed in Fig. 3.
Corresponding to the simulations in Sect. 4.2, we perform two cases of experiments: (1) a rotation around X axis of d = [5, 0, 0] T deg/s, and (2) a rotation around Z axis of d = [0, 0, 5] T deg/s. The desired trajectories of the angular velocity and  Figure 34 shows the control inputs of the system. In addition, the singularity index is indicated in Fig. 35. Those experimental results correspond to the simulation results in Figs. 24,6,7 and 27 in Sect. 4.2. Next, the experimental results of the Z-axis rotation are shown in Figs. 36,37,38 and 39. Similarly, the angular velocity, the vectorial part of the quaternion, the control inputs, and the singularity index are respectively shown in Figs. 36 37, 38 and The unmodeled gravity torque due to the misalignment between the center of rotation and the center of gravity of the testbed inevitably arises in the ground experiments unlike the simulations. In particular, the X-axis rotation is seriously affected by the gravity torque disturbance. Moreover, under this setting, since the singular direction of two of the four CMGs is aligned to that of the required torque, the X-axis rotation is much difficult than the Z-axis rotation. Nonetheless, the present MPC controller compensates the disturbance and successfully achieves attitude

Conclusions
Control Moment Gyros (CMGs) are widely used for precise attitude control. This paper describes the attitude control problem of a spacecraft, using a Model Predictive Control (MPC) systems. The features of the considered linear MPC are: (i) a virtual reference, to guarantee input constraints satisfaction, and (ii) an integrator state as a servo compensator, to reduce the steady-state error. Moreover, the real-time implementability is discussed, considering an optimized solver for embedded systems. Good performance are observed in simulations. However, the  proposed linear MPC control system is not able to handle with additive disturbances acting on the experimental testbed and it is only partially able to compensate this disturbance. In future works a control system able to handle and compensate the gravity gradient torque will be considered. Moreover, the gravity gradient disturbance will be also included in the model of the testbed itself.

Appendix: preliminaries on model predictive controller
Model Predictive Control is a control strategy that uses a dynamic model of the system, in order to make a prediction of the future behaviour of the system itself and consequently to "choose" the best control input. Since it is a modelbased approach, the definition of the mathematical model of the system is important to obtain a good control performance. From one hand, the model has to be descriptive enough to capture all the most important characteristics of the system, in order to make predictions as much as possible near to the real evolution of the process. On the other hand, if the model is too detailed, the complexity is increased, and the control system reaches a high computational cost, that can be processed just by expensive hardware devices Rawlings and Mayne (2009). At each sample time, the controller computes the sequence of control inputs for the actual and next sample times in open loop fashion, but just the first one of this sequence is applied to the plant. At the next sample time, the procedure is repeated from the beginning, without considering the other control input values obtained in the previous sample time. The number of control inputs, collected in the control sequence, depends on the capability of the controller itself to "look forward" in the future. This characteristic is defined by a parameter called Receding Horizon H s .
The following linear dynamic model of the system in state space formulation is considered: The first equation is the state equation, where x ∈ R n p is the system state, and u ∈ R m is the control input. The second equation is the output equation, and y ∈ R n y is the output vector.
The matrices of the state space system are: (1) the state matrix A ∈ R n p ×n p , (2) the control matrix B ∈ R n p ×m , (3) the output matrices C ∈ R n y ×n p and D ∈ R n y ×m . The notation x(k + 1|k) indicates the value of the state x at time k + 1 calculated at time k.
The prediction matrices S and T are defined using the defined model of the system, for the evaluation of the optimal sequence of control inputs (38) x(k + 1|k) = Ax(k|k) + Bu(k|k) (39) y(k|k) = Cx(k|k) + Du(k|k).
At each sample time a convex optimization problem is solved, where a cost function is minimized, to obtain the optimal sequence of the control inputs. The cost function at time k is defined as where Q,R and P are the weighting matrices, defined as follows The cost function can be rewritten as Matrix P is the terminal cost function, to ensure feasibility and stability to the MPC. Combining the prediction equation (40) and cost equation (42), we have

Rewriting the previous equation:
One of the main feature of the MPC controller is the definition of the input and state constraints, in the optimization problem. The input and state constraints are defined as follows: Finally we can define the constrained optimal control problem (Quadratic Programming) Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. (43)