C Rengifo

J.-R Chardonnet
email: jean-remy.chardonnet@ensam.eu

D Paillot
email: damien.paillot@u-bourgogne.fr

H Mohellebi
email: hakim.mohellebi@renault.com

A Kemeny

Solving the Constrained Problem in Model Predictive Control Based Motion Cueing Algorithm with a Neural Network Approach

Keywords: Optimization, recurrent neural networks, MPC motion cueing algorithm, real-time simulator, penalty method

Because of the critical timing requirement, one major issue regarding model predictive control-based motion cueing algorithms is the calculation of real-time optimal solutions. In this paper, a continuous-time recurrent neural network-based gradient method is applied to compute the optimal control action in real time for an MPCbased MCA. We demonstrate that by implementing a saturation function for the constraints in the decision variables and a regulation for the energy function in the network, a constrained optimization problem can be solved without using any penalty function. Simulation results are included to compare the proposed approach and substantiate the applicability of recurrent neural networks as a quadratic programming solver. A comparison with another QP solver shows that our method can find an optimal solution much faster and with the same precision.

Introduction

Driving simulators are necessary tools to evaluate human behaviors in driving situations, to validate benefits of advanced systems and autonomous driving and to simulate driving scenarios in safe and controlled environments. Motion commands such as accelerations and rotations cannot be sent directly from the model of the simulated vehicle to the platform, because of specific physical constraints and limits in the workspace.

For this reason, different control strategies called Motion Cueing Algorithms (MCA) are implemented, which aim to better reproduce all the signals coming from the simulated vehicle model while keeping the platform within its boundaries. There are mainly four motion cueing strategies, three of them being filter-based control strategies: classical, adaptative, optimal, and one based on prediction [Gar10]. Past works have shown that model predictive control (MPC)-based MCA [Dag04] is the most suitable strategy, as it uses optimization to track the reference, respects the constraints and maximizes the workspace. Moreover the MPC procedure provides easier tuning and is more robust [Neh06,Fan12,Beg12,Ven16].

In MPC-based MCA the system dynamics varies quickly and the numerical solution that is required in real time is often seen as an obstacle for normal solvers. This optimization problem solution in real time has been studied intensively [Goo06]. One promising approach is to use artificial recurrent neu-ral networks (RNN) as a method to compute dynamic optimization based on hardware implementation. The recurrent neural network capabilities such as high parallelism, adaptability, and hardware-based circuit implementation are recognized by researchers and are used in many different application areas including: control system design, associative memory, mathematical programming, signal and image processing.

In the optimization field, the first recurrent neural network model was proposed by Tank and Hopfield in their seminal work as an implementation on analog circuits [Tan86]. Since then many neural network models have been developed to solve different types of optimization problems. Some of them are summarized in Table 1.

Projection

Without penalty function.

[Xia05],[Che09],

[Liu15] Discrete time Good for hardware implementation.

[PI13]

Several studies have shown the neurodynamic optimization capability to solve quadratic problems relative to MPC strategy [Que93,Yun08,Yan12,Wan14].

To the best of our knowledge, none of these applications evaluates the applicability of the penalty method as a form of effective optimization for the implementation in an MPC loop.

In this paper, we propose to implement a continuoustime recurrent neural network based gradient method as a QP solver in the platform control of a dynamic driving simulator. A penalty method was proposed to solve the constrained optimization problem in real time. Compared to an active-set method [Goo06], we prove that, by using a saturation function for the decision variables and an appropriate regulation in the energy function, the penalty terms can be avoided and then the optimal solution gives a good performance in a tracking reference task in the MPC.

MPC-based MCA

Problem formulation

The dynamic driving simulator considered in this study is composed of an hexapod mounted on a XYtable . Since it has a dynamic that is predefined by the manufacturer, we consider the double integrator as a perfect model system. The idea with this model is to be able to control for each degree of freedom the longitudinal and rotational accelerations of the platform.

The linear time invariant system is described by the discrete-time state-space model:

x m (k + 1) = A m x m (k) + B m u(k) y(k) = C m x m (k) (1)
where x m (k) ∈ R s is the state vector: position, velocity, u(k) ∈ R l is the input vector: the reference acceleration, y(k) ∈ R t is the output vector: the simulator acceleration, A m ∈ R s×s and B m ∈ R s×l . The system is subject respectively to the following boundaries conditions:

x min ≤ x(k) ≤ x max y min ≤ y(k) ≤ y max ∆u min ≤ ∆u(k) ≤ ∆u max (2)
In a general tracking situation, the regulation error and the control action cannot be equal to zero at the same moment, thus we augment the model of the system with an embedded integrator to ensure offset free tracking in the control law [Ros03]. The input is the control increment ∆u(k) rather than the command u(k) itself:

x(k + 1) = Ax(t) + B∆u(k) y(k) = Cx(k) (3)
where,

x(k+1)

∆x m (k + 1) y(k + 1) = A A m 0 T m C m A m I m x(k) ∆x m (k) y(k) + B B m C m B m ∆u(k) y(k) = C [ 0 m I m ] ∆x m (k) y(k) (4)
The vectors of the predicted output Y and states X and the future control ∆U are created recursively and defined as:

X = [ x(k + 1) x(k + 2) • • • x (k + N p ) ] T Y = [ y(k + 1) y(k + 2) • • • y (k + N p ) ] T ∆U = [ ∆u(k) ∆u(k + 1) • • • ∆u (k + N u -1) ] T ( 5
)
where N p denotes the predictive horizon (1 ≤ N p ), N u denotes the control horizon (0 < N u ≤ N p ). In a more compact form:

Y = F x(k) + G∆U X = F x x(k) + G x ∆U (6)
where,

G x =       B 0 . . . 0 AB B . . . 0 A 2 B AB . . . 0 . . . . . . . . . . . . A Np-1 B A Np-2 B . . . A Np-Nu B       F x =     A A 2 . . . A Np     (7) F = CF x and G = CG x (8)
In the MPC also known as the receding horizon control, an open-loop control sequence is computed at each sampling time by optimizing an objective function over a finite predictive horizon N p , within the system constraints (Eq. 2). The first control action is implemented making a closed-loop system. The process is repeated at the next time step, yielding a receding horizon control strategy [Bor07]. The cost function to minimize is of the form:

J(k) = Np j=N1 δ(j) r (k + j|k) -y(k + j) 2 2 + Nu j=1 λ(j) ∆u(k + j -1) 2 2 + Np j=N1 q(j) x (k + j|k) 2 2 (9)
where J(k) is a scalar objective, subject to the constraints defined in Eq. 2, r (k + j|k) denotes the reference trajectory at sampling instant k, y(k + j) is the predicted output, δ(j), λ(j) and q(j) are weighting parameters for the tracking error, the control rate and the states respectively.

The prediction reference r(k + j) remains constant over the prediction horizon N p , overcoming the difficulty in predicting the driver's behavior.

In the matrix form, we get

J = (Y -R s ) T Q δ (Y -R s ) + ∆U T Q λ ∆U + X T Q q X where R T s = r(k) × [1 1 ... 1] 1×Np Replacing the predicted output Y = F x (k i ) + G∆U and the predicted states X = F x x (k i ) + G x ∆U in
J and eliminating all the terms that do not depend on the vector decision variable ∆U , the cost function obtained is

J = ∆U T G T Q δ G + Q λ + G x T Q q G x ∆U +2 (F x (k) -R s ) T Q δ G + (F x x (k)) T Q q G x ∆U
Finally, the QP problem with the constraints remains:

minimize J = 1 2 ∆U T H∆U + f T ∆U subject to A c ∆U ≤ b (10)
where,

∆U ∈ R n , f ∈ R n , b ∈ R m , A c ∈ R m×n
The matrix H is a symmetric positive definite n × n, i.e., all eigenvalues of H are strictly positive and then, the cost function J has a unique global minimum. The matrix equation A c ∆U ≤ b contains all the linear inequality constraints of Eq. 2. This quadratic problem has to be solved at every sampling time to find the control sequence ∆U and apply the first element of this sequence to the real system as a receding control strategy.

Feasibility of the solution

Many past studies have addressed the feasibility issue at each time step in the optimization problem due to the constraints set on the output and the state constraints. Different solutions exist to address this issue: making the horizon control infinite, setting a terminal state equal to zero, using a discrete linear quadratic regulator terminal condition. The first leads to an optimization problem with infinite constraints that is impossible to solve, the second can considerably reduce the feasible QP solution and might perturb the input trajectory in short horizons. The last one is difficult to apply to a real-time dynamic simulation system [Fan17].

Getting a feasible solution at every sample time requires an extremely high prediction horizon, which is not viable to find the solution in real time. Also previously applied methods are not suitable for MPCbased MCA. To overcome this issue, Fang et al. proposed a different condition [Fan14]. Their proposed approach can verify a solution along the prediction horizon using a braking law once the platform approaches its limits. This law can be compared to an adaptive filter that allows the platform to return to its neutral position while respecting the limits in position, velocity and acceleration. This condition is summarized as:

p min ≥ p k + c v T x + c a T 2 2 ≤ p max (11)
where the coefficients c v , c a and T are tuning parameters that prevent the platform from exceeding its limits. In Table 2 we expose the performance of the dynamic platform we consider here along the degrees of freedom of interest, here the longitudinal and lateral rail-axis.

Recurrent Neural Network model design

The idea of using optimization-based neural networks is to compute a QP optimal solution explicitly and to apply the control action online. The network used in this work is a continuous feedback network. Its dynamics can be depicted by a set of continuous dynamic systems as follows:

ẋ = w (x(t)) (12)
For a better understanding, we express the decision vector ∆U in Eq. 10 in terms of the states x in the network represented by Eq. 12.

minimize

f (x) = 1 2 x T Hx + f T x subject to g i (x) = A c x -b ≤ 0, i = 1, 2...m. ( 13
)
where f (x) is a strongly convex objective function since H is positive definite matrix, g(x) ∈ R m is the inequality constraints. The network represented by Eq. 12 needs a corresponding Lyapunov function E (x(t)) that defines the stability of the dynamic system trajectory thus, as t increases, the value of E (x(t)) decreases. It has been proven [Zha13] that the objective function of the optimization problem in Eq. 13 is equivalent to the Lyapunov function of the neural network.

The differential equations in Eq. 12 are solved simultaneously by an associated circuit consisting of highly interconnected processing units called neurons. For some initial state, the network will converge to an equilibrium state that will coincide with the optimal solution to the original problem (Eq. 13).

To solve the QP problem in Eq. 13 and to take into account the system constraints, the next step is to transform the constraint problem into an unconstrained problem to be able to apply a steepest descend method based on the gradient of the function. By using the penalty method, the energy function is adapted to integrate inequality constraints into the original function.

E(x, k) = f (x) + k m i=1 g i + (x) 2 g i + (x) = max {0, g i (x)} (14)
where k > 0 is the constant penalty parameter and g i + (x) is a continuous non-negative penalty function which equals zero at a point if only all constraints are satisfied at that point. The augmented energy function remains a quadratic objective function in which each unsatisfied constraint influences the state x by assigning a penalty multiplied by k. This parameter influences the accuracy of the solution; therefore, it must be large enough to provide a feasible solution. Ideally, when the value of k holds at infinity, the optimal solution of E(x, k) converges to a solution of the original problem with constraints (Eq. 13). In practice this remains impossible to realize since k cannot be set to infinite, therefore, with this method, the solutions are an approximation to the optimization problem (Eq. 13). If k is small, a greater violation of the restrictions will occur since the penalty is not enough to guarantee a feasible solution, on the other hand with a very large value of k the energy function will be computationally ill-conditioned and become a stiff system.

According to the literature [Coc93], it is desirable to create an activation function S i for the penalty function. That is, every time a constraint is active, the penalty function is different to zero, otherwise it's zero and the solution does not take into account the constraints and it's solve as an unconstrained optimization problem:

E(x, k) = f (x) + k 2 m i=1 S i [g i (x)] 2 (15)
In order to find the minimal value of the energy function E(x, k), the unsconstrained optimization problem is transformed into an associated system of firstorder ordinary differential equations by using the gradient descent method [Coc93]:

dx dt = -µ∇ x E(x, k) (16)
The state trajectories x(t) moves in the search direction of the negative of the gradient and will converge to a global minimum (equilibrium point in Eq. 16) for E(x, k), since E(x, k) is a Lyapunov function. The dynamic equation at neuron j can be written

dx j dt = -µ ∂f (x) ∂x j + k m i=1 S i g i (x) ∂g i (x) ∂x i , (j = 1, 2, ..., N u ), (i = 1, 2, ..., m), k > 0, µ > 0 (17) µ = diag (µ 1 , µ 2 , . . . , µ Nu ) , µ j = 1 τ , S i = 0, if g i (x) ≤ 0 1, otherwise
where µ is a learning parameter, τ is the integration time constant of the integrators, m is the number of inequalities constraints. Equations 16 and 17 represent a recurrent neural network as a gradient dynamic system. Replacing the functions by their respective indexes:

dx j dt = -µ j f j + Nu i=1 h ji x i + k m i=1 S i a ij g i (x) (18)
For convenience, the adaptive gains µ j are equal.

The selection of the appropriate parameter µ is crucial to ensure the stability of the system and the convergence speed to an equilibrium point. Considering the following aspects, the choice of µ is based on a trial and error basis: higher learning rates will decay the loss faster, but if it is too large, the algorithm will become unstable and steepest descent will not converge.

Implementation

To implement the proposed RNN we used Simulink-Matlab. As shown in the previous section, the recurrent neural network is seen as a set of differential equations ( 16) which represent all the state trajectories x of the network. A simplification of the general topology of the RNN can be seen in Fig. 1. The architecture consists in a continuous-time integrator that can be viewed as a special time-depending neuron; dynamic building blocks seen as multipliers or adders, function generators for the realization of the gradient of the energy and saturated transfer functions.

Figure 1: Scheme of the proposed one-layer recurrent neural network

In order to generate code and make a real-time based Simulink application, it is required to use a fixed-step integration algorithm to solve the system in Eq. 16. We simulate the RNN as a optimizer in the control loop as shown in Fig. 5, employing a penalty function method (Eq. 17) for inequality constraints. The method shows a convergence for all the state's trajectories x towards an optimal point since the solution tracking reference is the desired one Fig. 7. Nevertheless, when implementing the available fixed-step time solvers, the method did not work completely since the convergence of Eq. 17 was not ensured for all sample times.

Aiming to find the minimum step that the RNN needs to find a solution, i.e., an equilibrium point for every state x at every sample time, we analyze the step used by a variable-step solver. Fig. 2 shows that the minimum step size to capture the dynamics accurately is very small (∼ 10 -8 ) and cannot be applicable to a real time application. As seen in Fig. 2, the differential equation system is stiff, i.e., the system has both slowly and quickly varying continuous dynamics. Therefore it is mandatory to take small time steps in the numerical method to obtain satisfactory results. One of the reasons why the differential equations are stiff is that the system is ill-conditioned:

cond(H) >> 1 (19)
This condition means that the ratio of the largest eigenvalue to the smallest one of matrix H is very high. In this case, the gradient might not choose the best way to descend to the minimum and the solution to the QP problem may become unreliable. Fig. 3 shows the condition number of the Hessian matrix H as a function of the prediction horizon N p . However, the energy function does not depend only on the original cost function J (Eq. 13), but also on the penalty function (Eq. 14) that is updated at each sampling step according to the system's constraints g(x) and the states, thus making difficult the preconditioning of the system based solely on the Hessian matrix H of f (x).

Accordingly, to alleviate ill-conditioning and the stiffness of the differential equations, a regulation strategy is implemented. In Fig. 1 we can see a gain g that represents the regularization parameter that helps handling oscillations for each state trajectory x. Additionally, to ensure that all the constraints in Eq.13 are satisfied, and based on the convergence analysis in [Cos08], we convert the constraints for the system to boundary constraints for the control signal (∆U = x for the RNN). We choose a saturation linear function depicted in Fig. 4 that ensures the control does not exceed its boundaries:

x j = f (a j ) = l min , a j ≤ l min a j , a j ∈ [l min , l max ] l max , a j ≥ l max (20) j = 1, 2, ..., N u
This function will saturate each state trajectory x to an upper l max or lower l min limit based on the constraints and the actual dynamic states of the platform at each simple time. The RNN will find at a fixed-step time a control sequence ∆U = x for which the gradient of the cost function is zero. When the dynamic platform approaches the limits that are governed by the restrictions on the system, the saturation function will prevent the gradient from finding an unstable solution and will limit the value to the maximum or the minimum gradient as shown in Eq. 21.

x j = l max , ∂E(x) ∂x j ≤ 0 l min < x j < l max ∂E(x) ∂x j = 0 x j = l min , ∂E(x) ∂x j ≥ 0. (21)
The proposed dynamic system describing this behavior is detailed in the following state equations:

dx dt = -µ (∇ x E(x) + g.a) ( 22
)
Since the Hessian matrix H of the function f (x) is positive-definite, each equilibrium point of the system (Eq. 22) corresponds to an optimal global solution of (Eq. 13).

It is important to mention that the recurrent neural network used in this work (Fig. 1) has N u neuronsbased continuous time integrators that considerably improve the computational cost and allow us easy implementation. Fig. 5 shows the block diagram of the proposed control scheme.

Control Scheme

Based on the recurrent neural network design described above and making reference to the MPCbased MCA control strategy (Fig. 5), the control scheme can be summarized as follows:

1. Define the prediction horizon N p , the control horizon N u , the sample time, the weighting parameters δ(j), λ(j) and q(j) in Eq. 9. 2. Transform Eq. 9 in a QP problem (Eq. 13) and update the matrices H, f, A, b with the actual states and constraints. 3. Solve the quadratic programming problem using the recurrent neural network approach by calculating the optimal control actions ∆U . 4. Calculate the optimal control input u(k) = ∆u + u(k -1). 5. Apply only the first element in u(k) sequence and then apply it as an input to the system.

Simulation and results

Using the SCANeR Studio simulation software † , a driving simulation scenario was generated to take the vehicle information as the acceleration signal. The idea is to analyze the capacity of the simulator response to respect the limits and tracking the signal using the MPC-based MCA method with an optimization based on neural networks. The scenario consists in a simple maneuver remaining at an initial position for four seconds then following a car with a speed of 30 m/s in a straight line. The scenario lasts only 18 seconds which is enough to evaluate the behavior of the system as can be seen in the reference signal in Fig. 7. A fixed-time implicit solver was implemented, since it provides greater stability for oscillatory behaviors, avoiding unnecessary computational cost. The global sampling time of the system was set to 1 ms, and a minimum of two iterations are required for the network to converge to an optimal solution. The platform system is sampled with a 10 ms step that is the minimum required for a real-time application (control frequency at 100 Hz).

For the QP problem (Eq. 10), the tuning parameters were the followings: d = 1, λ = 0.1, q(pos) = 1, q(vel) = 0.1, q(acc) = 0, N u = 3, N p = 100. Optimization is implemented on each independent degree of freedom of the platform system, but here we only show the simulation of the most relevant axis: the longitudinal axis. First, we conducted a comparison to validate the accuracy of the optimal value at each time step. A quadprog-Matlab function is used as the comparison solver. This function can apply tree different algorithms: trust region reflective, interior point convex or active set. We decided to implement an active set method to find the QP optimal solution at each sample time. Using the previous information, we found that the error between the solution of that solver and the proposed approach based on RNN is almost null. In Fig. 6, we can see that the platform acceleration given by the control input of both curves overlap, proving the applicability of the RNN to solve the MPCbased MCA optimization problem. † http://www.oktal.fr/en/automotive/range-of-simulators/software Secondly, we analyzed the average time the quadprog takes to solve the optimization problem. This analysis was based using the model profiler application available in Matlab-Simulink. We found that for our problem with 3 states the quadprog function takes at least 3 ms to find a solution. In the simulink profile report can be seen the time taken by each block in solving the specific operation. For the optimization subsystem (Fig. 1) the total time for one sample time is 0.235ms. As stated above, a minimum of 3 iterations is required to find a optimal solution, then the time required to solve the Qp problem is ∼ 0.7ms. This comparison shows the potential of the proposal approach to solve problems with a greater number of states and solve the MPC-MCA obtimization problem in real time respecting the plateform constraints. It is observed that the platform tries to follow the reference of the acceleration until it approaches the limits. Then, the restrictions are applied so that the optimization can find a solution at each sampling time.

The way in which the acceleration behaves when it approaches the limits depends on the parameters mentioned in the feasibility section.

Comparing with previous recurrent neural networks based on penalty functions (Tab. 1) for solving a quadratic optimization problem with bounds constraints, our method shows to have lower computa-tional complexity since the penalty function can be replaced by using a saturated activation function for the bounds constraints and an appropriate regularization parameter. Consequently, the number of neurons and logical operations used to solve the QP problem can be reduced.

Conclusions and future work

The resulting set of differential equations requires less components with respect to existing neural networks, for the same class of optimization problems, and could converge to the optimal solution in a finite time. The simulation results confirm the effectiveness of our proposed method and better performance can be obtained with an appropriate learning rate tuning. Our method provides a good tracking performance with fast responses without exceeding the system limits in position, velocity and acceleration.

The neural networks parallel capability, gives a great opportunity for the implementation of this technique for a more complex system. Using a multi-core implementation, we can eliminate the adaptive filter condition of the constraints with a prediction horizon long enough to ensure optimization feasibility, avoiding the tuning task and making the algorithm much more robust.

As future works, we aim at tuning the learning rate so that the system can find a solution as quickly as possible within the required time and with an increased number of states. It is also intended to integrate the different axes of the platform, e.g., combine optimization between the longitudinal axis and the pitch angle. Future development will include a validation test of the algorithm on the real setup.

Figure 2 :

 2 Figure 2: Solver step time for a stiff model

Figure 3 :

 3 Figure 3: Condition number for the Hessian matrix of the system

Figure 4 :

 4 Figure 4: Saturation dynamic function f (a j )

Figure 5 :

 5 Figure 5: Block diagram of our proposed control scheme

Figure 6 :

 6 Figure 6: Comparison between quadprog and the proposed RNN

Figure 7 :

 7 Figure 7: Longitudinal acceleration rendering using rails only

Table 1 : Different RNN optimization approaches

 1

	Recurrent NN	Characteristics	References
	Penalty function	Easy	implementa-	[Maa92],[Mla00],
		tion, only approximate	[Che09]
		optimal solutions
	Lagrangian	Manage equality and in-	[Zha92],[Bou93],
		equality constraints.	[Cos08]
	Dual	Overcome	penalty-	[Xia96],[Shu06],
		method Lower compu-	[Hu08]
		tational complexity.

Table 2 : Performance along the longitudinal and lateral rails of Renault ULTIMATE simulator

 2

	Rail	Position	Velocity	Acceleration
	X	±2.6m	±2m/s	±5m/s 2
	Y	±2.6m	±3m/s	±5m/s 2