A * Conjugate of A ∈ C A T Transpose of A ∈ R A ≺ ( )0
Matrix A is symmetric and negative (semi)denite A ( )0

Matrix A is symmetric and positive (semi)denite Co(X)

Convex hull of set X A = A T Matrix A is real symmetric He(A) = A T + A

Introduction and problem statement of the thesis

In the recent years, study on advanced control methods for automotive systems has gained a huge momentum and most of the automotive industries has embarked in research and development and implementation of better control algorithms for improved passenger comfort and safety and as well as vehicle's performance in terms of energy eciency, pollution reduction etc. This sudden surge in interest is partly in response to the advent of autonomous vehicles and it has also been the key driving factor for automotive companies to strive for better and optimal performance. Conditioned upon the aforementioned requirements both objectively and subjectively, optimal control methods fares much better than other control methods due to the systematic approach embodied in its design to tackle the above issues. In order to implement optimal control in practice, model predictive control (MPC) formulation of the optimal control problem (OCP) provides the necessary feedback control framework to work seamlessly in real-time and real-world. MPC is one of the most ecient and powerful control methodologies and over the last few years, MPC has become commonplace both in industry and academia due to its performance and optimality. Initially, MPC (known as 2 Thesis framework and contribution dynamic matrix control) was restricted to chemical engineering with application for petrochemical industries [START_REF] Cutler | Dynamic matrix control?? A computer control algorithm[END_REF], however comprehending it's potential and performance benets, the method has gradually pervaded into other streams of engineering such as automotive, aerospace, biomedical, etc. Concomitantly, this has attracted several researchers from dierent engineering domains and this has positively ensued in multitude of its variants, methods, techniques and theory. Hitherto, some of the well-known extensions of the MPC [START_REF] Rakovi¢ | rndook of model preditive ontrol[END_REF] include (to name a few) explicit MPC [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF] (EMPC), nonlinear MPC [START_REF] Rawlings | Nonlinear model predictive control: A tutorial and survey[END_REF] (NMPC), Linear Parameter Varying MPC (LPV-MPC) [START_REF] Morato | Model predictive control design for linear parameter varying systems: A survey[END_REF], stochastic MPC [START_REF] Mesbah | Stochastic model predictive control: An overview and perspectives for future research[END_REF]] (SMPC), fault-tolerant MPC [START_REF] Nguyen | Fault Tolerant Predictive Control for Multi-Agent dynamical: formation reconguration using set-theoretic approach[END_REF]] (FTMPC), tube based MPC [START_REF] Mayne | Robust model predictive control of constrained linear systems with bounded disturbances[END_REF] (TMPC), learning based MPC [START_REF] Aswani | Provably safe and robust learning-based model predictive control[END_REF] (LMPC), economic MPC [Ellis, Durand, and Christodes 2014] (eMPC), adaptive MPC [START_REF] Bujarbaruah | Adaptive MPC for iterative tasks[END_REF]] (AMPC) etc. and more exciting extensions to be engendered in years to come and this trend seems to be growing unabatedly. The main reason for this popularity and rapid adoption stems from the idea of receding horizon control, the fact that an online optimization problem (in case of EMPC, an oine optimization problem is precomputed) is solved at every sampling period to obtain the optimal control inputs for the current state of the system. This provides the necessary leeway for the control engineer either to statically or dynamically include the required objectives and system constraints into the optimization problem for control.

With the given advantages of the MPC controller, in today's world, MPC is one of the highly coveted and sought method by the automotive industry due its ability to explicitly shape the performance and safety requirements by means of objective and constraint functions.

However, the major hindrance for the method from being pervasive is the need for high computational time for solving the online optimization problem. Over the last few decades with the development of advanced and sophisticated embedded processors, the gap between theory and practice for MPC controller has abridged. At the same time the MPC methods, techniques and its variants have evolved substantially in large proportions and number of research papers on these topics stands as a testimony for this fact [START_REF] Rakovi¢ | rndook of model preditive ontrol[END_REF].

Thus, it can be scarcely denied that with increased complexity of methods, a need for increased computational resource is inevitable to meet the real-time (RT) requirements, especially for fast sampled systems. Under this premise, the multicore hardware architecture of Graphic Processing Unit (GPU) displays a strong potential to provide a boost for unison of theory and practice for several state of the art MPC methods. GPU at its outset was primarily developed for accelerating the 3D graphics rendering pipeline for game engines and other multimedia applications. However in turn of events, due to its unique parallelism encompassed in its hardware (H/W) and software (S/W) architecture, it gained widespread attention in the scientic community and revealed a big deal in solving humongous simulations at faster time scales. In automotive domain, in the light of demanding computational needs, GPU based parallel computing has proven its mettle in solving RT path planning, navigation and control problems for autonomous vehicles [START_REF] Bojarski | End to end learning for self-driving cars[END_REF]].

The work proposed in this thesis combines the utility of MPC methodology as a control S/W paradigm and optionally, the H/W architecture of multi-core processors by the likes of GPUs and the proposed control scheme is termed Parameterized NMPC (pNMPC). In due course of the thesis work, a code generation S/W tool was developed and implemented in C++ environment with interface to MATLAB/Simulink for embedded implementation of the proposed control scheme for automotive systems. There are several automotive modules where control system play an important role in ensuring comfort and safety such as Electronic Stability Control system (ESC), Anti-lock Braking Systems (ABS), Automatic Emergency Braking (AEB), Collision Avoidance Systems (CAS) etc. [START_REF] Rajamani | ehile dynmis nd ontrol[END_REF]]. In this thesis the subject of focus is considered around modeling and control of suspension systems. The thesis was conducted under the purview of the EMPHYSIS project and the outcome of this thesis aligns with the goals of the project. The primary goal of the EMPHYSIS project is to utilize model based methods for design of advanced automotive control systems. The nal outcome of the project is to design and develop an automotive standard known as the Embedded Functional Mockup Interface (eFMI), which is used to deploy model based controllers onto the embedded automotive harware or Electronic Control Units (ECUs). GIPSA-lab is one of the academic partners of the project and the bestowed role is project demonstrator. Thus, GIPSA-lab is tasked with implementation of eFMU based controller for control of semi-active suspension system for the INOVE test platform and/or Hardware In the Loop (HiL) implementation in dSPACE MicroAutoBoxII (MABXII). This work is carried out in concerted eort along with its project partners from the French consortium which includes Siemens, SOBEN and Renault.

Vehicle suspension systems play a critical role in guaranteeing safety and comfort for the onboard passengers. There exists plethora of suspension systems depending upon the mode of operation and its technology, however under a bird's eye view, the entire spectrum can be briey classied into passive, semi-active and active suspension systems. Amongst the three classes, semi-active suspension systems are quite popular in the automotive industry due to multitude of reasons such as negligible power demand, safety, low cost and weight and signicant impact on vehicle performance [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]. Some of the prominent semiactive suspension technologies are a) Electro-Hydraulic (EH), b) Electro-Rheological (ER) and c) Magneto-Rheological (MH) based system. In this thesis, the main theme is considered around modeling and control of ER semi-active suspension system for the INOVE test platform [Vivas-Lopez et al. 2014a]. Thus, the aim of this thesis is to utilize the proposed pNMPC scheme for control of vehicle vertical dynamics by virtue of suspension systems.

In summary, the main contributions of this thesis are

To propose the pNMPC scheme for control of semi-active suspension system for control of vertical dynamics of vehicle. The method was experimentally tested and valided on the INOVE test platform for a quarter car system and also tested with Hardware In the Loop (HiL) simulations on dSPACE MicroAutoBox II (MABXII).

The thesis also proposes the plausibility of using multi-core processors by the likes of GPUs for solving the complex simulation based pNMPC scheme for fast sampled systems. As a proof of concept, this was tested in simulation for control of semi-active suspension system for a half car vertical dynamics model with the parallelized pNMPC scheme with INOVE parameters. The method was also augmented to incorporate the Thesis framework and contribution stochastic road prole model and this control scheme is termed as scenario stochastic pNMPC (SS-pNMPC). pNMPC -A code generation S/W tool was developed for solving a derivative free pN-MPC problem with the capability of generating both CPU and GPU codes. The code generation S/W was tested for a couple of examples for both CPU and GPU versions and also, tests (CPU version) on dSPACE MABXII were conducted for RT verication and validation of proposed code generation tool for a quarter car model with INOVE parameters.

Structure of the thesis

The organization of the thesis is presented with the following major parts:

Part I: Thesis and theoretical background Part II: pNMPC with RT applications for control of semi-active suspension system Part III: pNMPC -A code generation software tool for implementation of derivative free pNMPC scheme for embedded control systems

The rst part is dedicated to discuss the fundamentals and provide an outlook on the thesis and theoretical background.

Chapter 1 provides a general overview on automotive suspension systems and its classications, along with a brief overview on the mode of operation for each type. The vehicle vertical dynamics modeling is covered for quarter and half car models, which appear recurrently throughout the thesis. The performance and objective requirements of vehicle suspension systems are succinctly explained. The INOVE test platform's system congurations are expounded which includes the physical plant, software (S/W) and hardware (H/W) components. Finally, the chapter is concluded with a simplied explanation of the Functional Mockup Interface (FMI) standard and a brief overview on the EMPHYSIS project and its outcome -embedded FMI (eFMI) standard.

Chapter 2 provides the theoretical background for the thesis. The rst part of the chapter provides a quick overview on nonlinear optimization problem, its classications, sub-branches and followed by a literature survey on several o the self nonlinear optimization solvers. The second part discusses the fundamentals of NMPC methods and a brief literature survey on several MPC/NMPC toolboxes and nally, the third part is concluded with a gentle introduction to GPGPU computing (General Purpose computing on Graphic Processing Unit) with introduction to parallel programming with CUDA, CUDA memory model, atomics and other CUDA libraries.

Structure of the thesis 5

The second part discusses about the pNMPC scheme, its application and RT implementation. Also, the implementation of GPU based parallelized pNMPC scheme is explained in detail in Chapter 4.

Chapter 3 rst part introduces the concept of parameterized NMPC scheme and its application for control of semi-active suspension system for the INOVE test platform.

The pNMPC controller's design requirements for both objective and constraints for the suspension control problem are explained in detail. The second part explains the experimental implementation of the pNMPC controller for control of semi-active suspension system. The experimental study conducted to model the ER semi-active damper characteristics and parameter estimation are discussed in detail. From the obtained model parameters, the pNMPC controller is experimentally validated on the INOVE platform and also with HiL simulations conducted on dSPACE MABXII embedded target. The core results of this chapter are based on the following papers [Rathai et al. 2018], [Rathai, Alamir, and Sename 2019].

Chapter 4 discusses the GPU simulation implementation of the pNMPC controller for control of semi-active suspensions system in MATLAB/Simulink for the INOVE test platform model and parameters. In the rst part, the parallelized pNMPC scheme is introduced and in the second part, a stochastic version of the parallelized pNMPC is introducted and also, the method was tested in RT on several NVIDIA embedded boards to verify and validate the feasibility of the proposed scheme. The core results of this chapter are based on the following paper [Rathai, Sename, and Alamir 2019].

The third part is completely dedicated for the pNMPC code generation S/W tool, its working, features, implementation and along with few examples on CPU (PC), GPU and embeddeded implementation in dSPACE MABXII.

Chapter 5 introduces the pNMPC code generation S/W tool for implementation of a derivative free pNMPC scheme for embedded control systems. This chapter expounds the features of the developed code generation S/W and also, explains the code generation process in a bird's eye view. With a variety of examples, the proposed S/W was compared against ACADO toolkit and the simulation results are presented. Furthermore, using the proposed S/W, HiL tests were conducted on dSPACE MABXII for control of semi-active suspension system for a quarter car model with INOVE parameters and the simulation results are presented. Results for the GPU version of the code generation codes are also presented in this chapter and nally, the chapter is concluded with future works and conclusions. The pNMPC code generation S/W tool C/C++, MATLAB/Simulink codes are present in the following GitHub repository [START_REF] Rathai | pxwgX e gode qenertion ool por smplemenE ttion of pxwg gontroller por imedded gontrol ystems[END_REF]]. to minimize the roll angle of the vehicle (ride handling) [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]]. The physical realization of a suspension system is materialized with three main parts and an optional part which are a) the structure of the suspension system, which denes the suspension geometry, b) the spring element, which provides proportional and opposite force to the displacement of the suspension system and also acts as an energy storage element, c) the damper element aka shock absorber which provides proportional and opposite force to the velocity of the suspension system and also used to dissipate energy [START_REF] Goodarzi | Vehicle suspension system technology and design[END_REF], [START_REF] Gillespie | pundmentls of vehile dynmis[END_REF]], and d) the optional part, which involves a motion actuator which applies only for active suspension system, which is used to provide additional energy by exerting force on the system.

Part I Thesis and theoretical background

A typical anatomy of a suspension system is illustrated in Fig. 1.1. The suspension system is housed between the wheel (unsprung mass) and the chassis (sprung mass) of the vehicle and 1.2. Automotive suspension system 11 it supports the whole weight of the vehicle. The suspension system also play an important role in describing the lateral, longitudinal and vertical dynamics of the vehicle as the forces and moments along these respective axes are directly aected through the wheels and the suspension system links the wheels to the road, the chassis to the wheels, thereby directly transmits these eects to the chassis. These eects are well pronounced during cornering or evasive maneuvers where vehicle handling becomes an issue of concern.

On the basis of suspension mechanics, the elements that constitute the suspension system can be divided into two parts which are the dynamic elements and kinematic elements. In this thesis, the prime focus lies on the study and control of suspension systems based on the dynamic elements which includes the spring, shock absorbers etc. whereas on the other side of the spectrum i.e. the kinematic elements, which involves the suspension geometry, suspension mechanism etc. and the kinematic eects are not considered in this work. Furthermore, on the basis of technological and dynamic perspective, suspension systems can be broadly classied

Passive suspension system

The passive suspension system is one of the most basic and perhaps a ubiquitous system deployed in several commercial automobiles. Typically, a passive suspension system consists of a spring and a shock absorber and the system is void of any element of control for manip-ulating the suspension characteristics online [START_REF] Miller | Tuning passive, semi-active, and fully active suspension systems[END_REF]], [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]. The design parameters such as stiness coecient, damping coecient, suspension structure etc. are optimized and tuned during the research and development phase of the product and the design parameters are chosen based on several requirements such as road condition, rattle space requirements, maximum number of passengers etc. Typically, the spring element is constituted either by coil springs or gas springs which acts as the stiness component and the dissipative element is constituted by hydraulic shock absorbers which exhibit nonlinear characteristics such as hysteresis, friction etc. The pros of a passive suspension system are: low cost, low weight, free of any electronics or sensor components and the cons are: not suitable for varying road conditions where comfort and safety are of paramount importance and the suspension system is inexible to adaptation. The SER plot for passive system is pictorially illustrated by the red dashed line in Fig. 1.2 which lies in the rst and third quadrant in the SER plot, which typically has non-linear characteristics.

Active suspension system

In contrast to the passive suspension system, the active suspension system can generate, store and dissipate energy which is illustrated in Fig. 1.2, where the SER characteristics lies in all four quadrants [START_REF] Tora | The active suspension of a cab in a heavy machine[END_REF]]. Typically, an active suspension system consists of a spring, a shock absorber and an active actuator to impart energy to the system [START_REF] Fischer | Mechatronic semi-active and active vehicle suspensions[END_REF].

The actuator modulates the displacement between the chassis and wheels by providing input force which ought to be controlled [START_REF] Sun | edvned ontrol for vehile tive suspension systems[END_REF]. The pros of active suspension systems are: better ride quality, handling and safety and the cons are: exorbitant price, increased power consumption, heavy mechanical components and uncertain management of safety issues [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]].

Semi-active suspension system

Semi-active suspension system is partly an active suspension system, however, the system can not impart additional energy to the system and can only dissipate energy by manipulating the damper characteristics online. The region of operation is illustrated in Fig. 1.2 where the rst and third quadrant represents the semi-active damping region. However, the range of operation is much wider than passive suspension system and this provides the necessary latitude for controlling the system characteristics on a broader range. The pros of semi-active suspension system are: negligible power demand, safety characteristics, signicant impact on vehicle performance and low cost and weight of the system [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]] and the cons are: reduction in vehicle handling performance, ride quality when compared to active suspension system. Nevertheless, the system is adaptable to road conditions with better controllers. In this thesis, the focus is on control of semi-active suspension system and its application for control of vertical dynamics of vehicle. Also, this thesis follows the work of former PhD theses [START_REF] Poussot-Vassal | Commande robuste LPV multivariable de chassis automobile[END_REF], [START_REF] Do | LPV approach for the vehicles dynamics robust control: joint comfort and safety improvement[END_REF]] and [Nguyen 2016a] on control of semiactive suspension systems.

Vertical vehicle dynamics modeling 13

The semi-active suspension system could be classied based upon its technology and some of the most common variants are a) Electro-Hydraulic (EH), b) Magneto-Rheological (MR) and c) Electro-Rheological (ER) damping models. Despite the working mechanism of the semiactive suspension system for dierent technologies are distinct from each other, the underlying principle of operation, properties and characteristics are virtually the same for all semi-active damper classes. A brief insight into these working principle of these classes are provided below Electro-Hydraulic (EH) damper -The EH damper system consists of electronic valves, which are typically controlled by manipulating the solenoid spool valve and therefore, this operates the damper chambers for the uid to navigate [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]], [START_REF] Spelta | Design and applications of semi-active suspension control systems[END_REF]. By varying the uid levels, the damping coecient is controlled continously and linearly with the area of the opening valve [START_REF] Aubouet | Modélisation et commande de suspensions semi-actives SOBEN[END_REF]].

Magento-Rheological (MR) damper -The MR damper consists of Magneto-Rheological uid which is formed by suspending magnetic particles into oil. By varying the magnetic eld, the particles aligns itself to the eld and this in turn varies the viscocity of the uid [START_REF] Rabinow | The magnetic uid clutch[END_REF][START_REF] Lozoya-Santos | A toolkit for nonlinear model predictive control using gradient projection and code generation[END_REF]].

Electro-Rheological (ER) damper -The ER damper shares a similar principle of MR damper, however the damper uid is a mixture of oil and micron sized electric eld sensitive particles. Thus, by varying the electric eld, the damping co-ecient is changed and this principle is used for control of the damper [START_REF] Winslow | wethod nd mens for trnslting eletril impulses into mehnE il fore[END_REF]]. The semi-active dampers utilized for the INOVE test platform are ER dampers (for all four wheels of the vehicle) [START_REF] Nguyen | Semi-active suspension control problem: Some new results using an LPV/H∞ state feedback input constrained control[END_REF]] [Vivas-Lopez et al. 2014b].

Vertical vehicle dynamics modeling 1.3.1 Quarter car vehicle model

The quarter car model illustrated in Fig. 1.3 consists of two mass elements which are the sprung mass element (chassis) and the unsprung mass element (wheel). The vertical dynamics model for the system around equilibrium [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]] is expressed with

m s zs = -k s (z s -z us ) + u m us zus = k s (z s -z us ) -u -k t (z us -z r ) (1.1)
where, m s and m us are the sprung and unsprung masses respectively, k s and k us are the stiness coecients of the damper system and wheel respectively. z s , żs , z us and żus are the sprung mass position, velocity and unsprung mass position, velocity respectively and z r is the vertical road displacement. It is worth to notice that u in this model is the force exerted due to the semi-active suspension system. żdef = żus -żs and z def = z us -z s are the deection velocity and deection position between the sprung and unsprung mass. The dynamics equation (1.1) can be compactly expressed in state space form with

ẋ = A c x(t) + B c u(t) + B dist c d(t) (1.2)
A c =      0 1 0 0 -ks ms 0 ks ms 0 0 0 1 0 ks mus 0 -(ks+kt) mus 0      , B c =      0 1 ms 0 -1 mus      , B dist c =      0 0 0 kt mus      (1.3)

Half car vehicle model (Roll oriented model)

The half car vertical dynamics roll oriented model, illustrated in Fig. 1.4 is a 4 degrees of freedom (DOF) model [START_REF] Mahala | Mathematical models for designing vehicles for ride comfort[END_REF]Deb 2009] [Savaresi et al. 2010], which involves chassis dynamics (heave motion), roll dynamics and dynamics of the two unsprung masses (wheels). The model ought to be viewed as the vehicle being scrunched from the front and rear ends into a single block. Let the left and right corner of the vehicle be indexed with i ∈ {l, r} respectively. The 4 DOF mathematical model is expressed with the following equations

           m s zs = -i∈{l,r} F s,i I x θ = (l l F s,l -l r F s,r ) m us,l zus,l = (-F s,l + F t,l ) m us,r zus,r = (-F s,r + F t,r ) (1.4)
where, m s , m us,l , m us,r represent the chassis mass, unsprung masses for the left and right corners. I x represents the moment of inertia along the roll axis. l l and l r represents the length of the chassis from the left and right corners with respect to centre of gravity (COG). F s,i represents the chassis forces and F t,i represents the wheel forces ∀i ∈ {l, r} which are expressed

F s,i = -k s,i (z s,i -z us,i ) + u i F t,i = -k t,i (z us,i -z r,i ) (1.5)
where, k s,i and k t,i represents the stiness coecent of the semi-active suspension system and wheel respectively. z r,i and z us,i represents the vertical road displacement and unsprung mass position ∀i ∈ {l, r}. u i represents the actuation force ∀i ∈ {l, r} exerted due to the semi-active suspension system. z s,i , ∀i ∈ {l, r} represents the sprung mass displacement at each corner which are obtained from the following equations

z s,l = z s + l l sinθ z s,r = z s -l r sinθ (1.6) Let X = [z s , θ, z us,l , z us,r
, żs , θ, żus,l , żus,r ] denote the state vector (X ∈ R 8 ), U = [u l , u r ] denote the input vector (U ∈ R 2 ) and D = [z r,l , z r,r ] denote the disturbance vector (D ∈ R 2 ), then the half car model (1.4) can be compactly expressed using the nonlinear state space equation with

Ẋ(t) = f (X(t), U(t), D(t)) (1.7)

Performance and objective requirements for vehicle suspension control

Despite the fact that the suspension system plays a pivtol role in shaping the lateral, longitudinal and vertical dynamics of the vehicle, automotive engineers study and view the system Chapter 1. Thesis background predominantly through the lens of vertical dynamics as the eects of road roughness are very palpable on the vehicle body [START_REF] Goodarzi | Vehicle suspension system technology and design[END_REF]. The three important performance and objective requirements for suspension control are:

Ride comfort -Ride comfort objective is related to the comfort of the on-board passengers of the vehicle. The comfort objective can be quantied by vibration isolation due to road roughness. The spectrum of frequencies from 0 -20Hz is considered as the comfort zone and ideally, the vertical displacement of the chassis of the vehicle must be same of the road in low frequencies (< 1Hz) and attentuated for high frequencies (> 1Hz) [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]. This is achieved by minimizing the vertical acceleration or velocity or displacement of the chassis of the vehicle.

Road holding -Road holding objective requires the vehicle to remain in contact with the road which is highly sought for safety requirements of the vehicle. The cornering, traction and braking abilities depend on the lateral and longitudinal forces of acting on the tyre which is directly a function of the normal force [START_REF] Goodarzi | Vehicle suspension system technology and design[END_REF].

Thus, the primary goal of this objective is to minimize the uctuations in the normal force of the wheel and this is assured by minimizing the deection between the wheel and road position.

Ride handling -Ride handling is a general requirement which necessitates the vehicle to be stable in every maneuver and also an important condition for stability. The fundamental requirement of ride handling objective is to ensure that the vehicle behavior is predictable and this information is communicated to the driver. The suspension system play a crucial role for ride handling by minimizing the vehicle's roll and pitch motion, controls the wheel angles and decreases the lateral load transfer during cornering maneuvers [START_REF] Goodarzi | Vehicle suspension system technology and design[END_REF].

INOVE test platform (GIPSA lab)

The The physical plant is a 1:5 scaled baja styled racing car which refects a miniature version of a full vehicle which encompasses wheels, engine and braking system. The key component of the plant is the semi-active damping system.

The semi-active suspension system involves four ER dampers with a force range of ±50N . The damper is manipulated by varying the input voltage from 0V to 5kV which is in turn manipulated by varying the pulse width modulation (PWM) duty cycle (DC) signal through CarCon2 driver module. The frequency of PWM-DC signal is 25kHz.

The Remote Controlled (RC) car is mounted over four OMRON linear servomotor module over all four corners and this is controlled appropriately to generate user dened input road prole. The servomotors has a bandwidth of 0 -20Hz and a maximum linear velocity of 1.5 m/s.

Sensor description

To capture the complete vertical dynamical behavior of the vehicle, the plant is probed by means of several onboard sensors.

Chapter 1. Thesis background

To measure the vertical acceleration of the unsprung masses (wheels), four Texense 1-axis capacitive accelerometer are used.

The stroke deection (i.e. the displacement between the sprung mass and unsprung mass) and the road displacement (motor position) are measured using eight Gefran resistive linear displacement sensors for each corner of the platform.

To measure the unsprung mass displacement, four Micro-Epsilon draw-wire displacement sensors are used for all four corners of the platform.

To measure the acceleration and angular velocity of the chassis a SBG MEMS based Attitude and Heading Reference System (AHRS) unit is placed on the sprung mass of the plant. The AHRS unit measures 3 accelerations along lateral, longitudinal and vertical axes and 3 angular velocities along roll, pitch and yaw axes.

The force exerted by the ER damper is measured via four force sensors mounted on ER dampers across all four corners.

To measure the tire forces across all four corners, four force sensors are placed at the base of the 4-poster testbed. The functional mockup interface (FMI) is a tool independent standard used for model exchange and co-simulation of dynamical systems. The standard was the result of ITEA2 MODELISAR project, where FMI 1.0 was primarily developed for model exchange purpose [START_REF] Blochwitz | The functional mockup interface for tool independent exchange of simulation models[END_REF]. Ever since the release of its initial version, its popularity surged and several tool vendors adopted the standard and this paved way to FMI 2.0 which is used for both model exchange and as well as co-simulation purpose with improved and new features to ease the use and increase the performance for larger models [START_REF] Blockwitz | Functional mockup interface 2.0: The standard for tool independent exchange of simulation models[END_REF]]. An instance of FMI is called as Functional Mockup Unit (FMU), i.e. the FMI standard denes the blueprint and the FMU is a container which consists of all the source code les (C-code), meta-data model les (XML les), binaries (static/dynamic libraries) les which are compliant with FMI standard to represent the underlying dynamical model. As mentioned previously, the FMI standard consists of two parts: The model developed in a FMI compliant tool is exported to FMU, which encompasses the details of the dynamic system such as states, inputs, parameters etc. The simulation tool utilizes this information from the FMU and simulates it by means of its inbuilt solvers such as explicit/implicit integrators, xed/variable step solvers etc. It is also important to note that the quality of solution is contigent over the type of solver used 20

Chapter 1. Thesis background to simulate the FMUs and this in turn is provided by the simulation environment.

FMI for co-simulation -The general concept of co-simulation is to simulate multiple interdependent subsystems in the simulation environment. The FMU exported for cosimulation purpose is shipped with the model and as well as the solver to simulate the system independently. Typically, the simulation tool would consist of multiple FMU blocks which interact with each other through input/output (I/O) channels and the primary task of the simulation tool is to handle the communication trac between the I/O ports of the subsystem, synchronization of data transfer, signal extrapolation and error control. In short, the simulation tool micromanages the whole collection of FMU subsystems and simulates it holistically and seamlessly to reect the dynamics of the real-world system.

A schematic of the FMUs subsystem is shown in Fig. 1.7, where the block is assumed to be interacting with the simulation environment (For more details, refer [START_REF] Blochwitz | Functional mock-up interface for model exchange and cosimulation[END_REF]).

1.7 EMPHYSIS project and embedded FMI (eFMI) standard

Conclusions

In this chapter a brief overview on automotive suspension system is provided which involves the classication of suspension systems and more specically on classication of ER semiactive damper system based upon its technology. The two vertical vehicle dynamics modelsquarter and half car models are presented, which are pertinent to this thesis and recurs in the upcoming chapters. A qualitative explanation of the objective requirements for the suspension system is discussed and this notion is concretized in the following chapters. A quick overview of the INOVE experimental platform is provided and nally, the chapter concludes with a simple explanation on the FMI standard and also, the gravitas for the EMPHYSIS project as well its outcome i.e. the eFMI standard is expressed.

Introduction

This chapter is dedicated to provide the fundamentals on non-linear optimization problems, dierent methodologies for implementation of non-linear model predictive control (NMPC) schemes and nally, the chapter is concluded with a brief outlook on general purpose graphic processing unit (GPGPU) computing using CUDA programming framework.

Nonlinear Optimization

Introduction

A general formulation for an optimization problem is dened by the following form minimize

x f (x) subject to x ∈ X ⊆ X (2.1) 26
Chapter 2. Theoretical background where x is the vector of optimization or decision variables of interest and X is the admissible or feasible subset on the optimization domain (X) [START_REF] Borrelli | reditive ontrol for liner nd hyrid systems[END_REF]. The optimization problem is carried out with respect to an objective function f : X → R which is used to assign every vector x a cost value with f (x) ∈ R. Let the optimal cost value for the problem (2.1) be represented with f * and the optimizer or optimal solution with x * , then the following equation holds

f (x) ≥ f (x * ) = f * ∀x ∈ X , with x * ∈ X (2.2)
In simple words, equation (2.2) states that among all the possible feasible values of x, the solution x * provides the optimal solution for the objective function f (x). The optimal objective value f * is also dened as the greatest lower bound for the problem (2.1). If f * = -∞, then the problem is regarded as unbounded below and if f * = ∞, then the problem is regarded as infeasible and if X = X, i.e. the whole search space, then the problem is regarded as an unconstrained optimization problem. The problem of determining the existence of optimal solution is called as a feasibilty problem. The set of optimizers x * (the solution could either be singleton or set-valued) that optimizes problem (2.1), i.e. f * = f (x * ) is dened with the following form

argmin x∈X f (x) = {x ∈ X : f (x) = f * } (2.3)
The optimization problem dened in equation (2.1) is generic in nature, however depending upon the structure of the objective function f (x) and the search space domain X , there exists plethora of classications such as linear and non-linear, deterministic and stochastic, continuous and discrete, convex and non-convex, nite and innite dimensional and also, the combination of all the aforementioned divisions. In this thesis the prime focus lies on a specic class of optimization problem known as non-linear optimization or non-linear programs (NLP), which is dened by minimize

x∈R n f (x) subject to g(x) = 0, h(x) ≥ 0 (2.4)
where f : R n → R, g : R n → R p , h : R n → R q are the objective function, equality constraints and inequality constraints respectively and the involved functions are assumed to be continuous and smooth in nature. In case when the knowledge of these functions are available, one can make further assumptions such as the functions are once or twice dierentiable. In other cases, such as black box models which might involve computer codes, compiled binaries, eFMU/FMUs models or functions etc., only the function values could be queried and the derivative information may not be available for solving the pertinent NLP problem. Thus, it is important to address the nonlinear optimization problems for both these cases which appears in several real-world applications. Typically, all the methods utilized to solve the NLP problem (2.4) are iterative in nature and the solution converges either to the optimal or suboptimal values depending upon various factors such as accuracy, relative tolerance, number of iterations etc., which are the tuning parameters for the NLP solver.

Classication of nonlinear optimization problems

As mentioned previously, in this thesis the prime focus lies on nonlinear optimization problems with the assumption that the involved functions are continuous and smooth in nature. Under this condition, the whole domain of NLP could be broadly classied into convex and nonconvex problems. The key property of a convex problem is that the local optimizer is the global optimizer, thereby the obtained solution is the global optimimum for the optimization problem.

Convex optimization involves optimization of convex functions over convex sets [START_REF] Boyd | gonvex optimiztion[END_REF] only function values can be queried for a given argument list or input vector. Typically, it is always recommended to utilize the derivative information whenever it is available, however, in situations where the details of the model or the functions are concealed from the end user due to security reasons or privacy reasons or protection of intellectual property rights, then it becomes highly impractical and cumbersome to utilize derivative based optimization methods.

Thus, the natural recourse is to resort to derivative free methods for solving NLP problems.

A brief outlook into derivative based and derivative free methods are provided in the following sections. It is important to note that the primary intention of this chapter is only to provide an overview and not delve into the implementation details of any of the methods.

Nevertheless, appropriate references are cited to benet the reader to study more in-depth on these topics.

Derivative based methods

The state of the art derivative based methods can further be classied into the following types First order derivative (Jacobian/gradient) based methods

Second order derivative (Hessian) based methods

Quasi-second order derivative based methods First order derivative (Jacobian/gradient) based methods: Over the last few years, the popularity of gradient based methods has increased in several folds and the main reason for this sudden surge in interest is attributed to recent developments in the eld of machine learning (ML) [START_REF] Goodfellow | ivluting derivtivesX priniples nd tehE niques of lgorithmi di'erentition[END_REF]]. When dealing with suret of training data, the rst order derivative based methods provide faster and better solutions for large scale ML problems.

There exists variety of gradient based methods such as classic gradient descent, batch gradient descent, stochastic gradient descent, gradient descent with momentum, Nesterov accelerated gradient (NAG), adagrad, adadelta, RMSprop, adams, adamax etc. The survey paper [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF]] provides a quick glimpse into the several gradient based methods utilized for solving NLP problems in the ML community. However, the aforementioned gradient based methods are mostly designed for unconstrained NLP problems. In optimal control, one is obliged to include the constraints into the NLP problem formulation. Some of the constrained gradient based methods include proximal gradient methods [START_REF] Parikh | Proximal algorithms[END_REF], Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | histriuted optimiztion nd sttistil lerning vi the lternting diretion method of multipliers[END_REF], operator splitting methods [START_REF] Stathopoulos | Operator splitting methods in control[END_REF] etc. and also, several variants of these methods. These are commonly adopted for online optimization in optimal control community. The prime benet of rst order method is, it is not computationally taxing compared to other methods, however, the major downside is that the convergence rate is either superlinear or linear, which is slow when compared to second or quasi-second order methods.

There are several rst order derivative based solvers for solving the unconstrained NLP problems and these are predominantly utilized in the ML community [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF]]. However, for the constrained case and as well as in the context of optimal control and RT implementation, this is currently in its inchoate state and there are only a handful of promising solvers such as OSQP [START_REF] Stellato | OSQP: An operator splitting solver for quadratic programs[END_REF], FiOrDOs [Ullmann 2011a], TFOCS [Becker, Candes, and Grant 2011], QPgen [Giselsson 2018], GPAD [START_REF] Patrinos | An accelerated dual gradient-projection algorithm for linear model predictive control[END_REF], POGS [START_REF] Fougner | Parameter selection and preconditioning for a graph form solver[END_REF], SCS [O'donoghue et al. 2016], FORCES PRO [START_REF] Zanelli | FORCES NLP: an ecient implementation of interior-point methods for multistage nonlinear nonconvex programs[END_REF]. For more details into the inner workings of rst order methods, refer [START_REF] Beck | sntrodution to nonliner optimiztionX heoryD lgorithmsD nd ppliE tions with wevef[END_REF]].

Second order derivative (Hessian) based methods: The whole class of second order derivative methods relies on multiple variants of the classic Newton's method for solving the roots of non-linear equations. In short, the NLP formulation in (2.4) is converted to a set of nonlinear equations by means of KarushKuhnTucker (KKT) conditions and at every Newton iteration of the NLP problem, a linear algebra subroutine is solved and this is repeated until convergence to the optimal solution. Despite the simplicity of the above explanation, the implementation methods are quite involved and requires special techniques such as matrix pre-conditioning, determining the right step size, computational complexity etc. Programming (IPM-QP) methods. A detailed explanation on the aforementioned methods are discussed in [START_REF] Bazaraa | xonliner progrmE mingX theory nd lgorithms[END_REF]. Typically, second order methods possesses quadratic rate of convergence and the optimal solution is generally obtained within a couple of iterations. Despite the fast convergence of the second order methods, one of the main issues is the demand for high computational need to compute the Hessians and this doesn't scale well for large scale problems. However, in optimal control, the structure of the NLP problem is exploited due to the specialized functional form of the objective and constraints functions and the Hessians are pre-computed and either cached or stored oine (the Hessians are mostly sparse matrices). This provides a signicant boost for RT implementation of optimal control problems in embedded systems.

Quasi-second order derivative based methods: As mentioned previously, one of the major issues with second order methods is the high computational complexity involved in computing the Hessian matrices. The conventional methods used to compute the derivatives for Hessians are nite dierences method and algorithmic dierentiation. However, the nite dierences method is highly prone to numerical inaccuracies and algorithmic dierentiation may not be possible for implementation for large scale NLP problems. Quasi-second order methods set a middle ground by retaining the computational complexity of rst order methods and convergence property of second order methods [START_REF] Dennis | Quasi-Newton methods, motivation and theory[END_REF]. The quasi-second order methods approximates the Hessians using the zeroth order information and rst order derivative (Jacobians) information and these approximations are plugged into Chapter 2. Theoretical background the IPM or SQP solver. The approximation is improved at every Newton iteration and this procedure is repeated until convergence to optimal solution for the NLP problem. Some of the well known Hessian update rules under the quasi-second order methods are Broy-denFletcherGoldfarbShanno (BFGS) algorithm, Broyden's rank-1 update and family of its method, DavidonFletcherPowell (DFP) algorithm, Symmetric rank-one methods etc. A detailed exposition into these methods are provided in [START_REF] Nocedal | xumeril optimiztion[END_REF].

The list of commonly used second order and quasi-second order NLP solvers are listed in Table 2.1. For more details on the classication of the NLP solvers, refer [START_REF] Leyer | Nonlinear constrained optimization: methods and software[END_REF]. The references for the respective solvers are: CVXOPT [START_REF] Vandenberghe | The CVXOPT linear and quadratic cone program solvers[END_REF]],

IPOPT [Wächter and Biegler 2006a], KNITRO [START_REF] Byrd | KNITRO: An integrated package for nonlinear optimization[END_REF], LOQO [START_REF] Vanderbei | LOQO: An interior point code for quadratic programming[END_REF]], CONOPT [START_REF] Drud | CONOPTa large-scale GRG code[END_REF]], FilterSQP [START_REF] Fletcher | User manual for lterSQP[END_REF]], LINDO [START_REF] Lin | The global solver in the LINDO API[END_REF], LRAMBO [START_REF] Griewank | Maintaining factorized KKT systems subject to rank-one updates of Hessians and Jacobians[END_REF], NLPQLP [START_REF] Schittkowski | NLPQLP: A fortran implementation of a sequential quadratic programming algorithm with distributed and non-monotone line searchuser's guide[END_REF]], NPSOL [START_REF] Gill | ser9s guide for xyv @version RFHAX e portrn pkge for nonliner progrmmingF Tech[END_REF], SNOPT [START_REF] Gill | SNOPT: An SQP algorithm for large-scale constrained optimization[END_REF], SQPlab [START_REF] Gilbert | l!e wtl softwre for solving nonliner optimiztion proE lems nd optiml ontrol prolems[END_REF]], qpOASES [Ferreau et al. 2014a], FORCES PRO [START_REF] Domahidi | FORCES Professional. embotech GmbH[END_REF]. 2006], [START_REF] Beck | sntrodution to nonliner optimiztionX heoryD lgorithmsD nd ppliE tions with wevef[END_REF]], [START_REF] Larson | Derivative-free optimization methods[END_REF], [START_REF] Conn | sntrodution to derivtiveEfree optimiztion[END_REF]. A comparative study of several state of the art derivative free optimization methods are presented in [START_REF] Pham | Comparative study of derivative free optimization algorithms[END_REF], [START_REF] Moré | Benchmarking derivative-free optimization algorithms[END_REF], [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF]].

It is also important to note that the eld of derivative free optimzation is not as mature as derivative based methods and most of the research in the past have been conducted only for the unconstrained case. Only a few methods exists for constrained case and this eld is still being researched for new directions. Table 2.2 lists the commonly used derivative free optimization solvers, for more details refer [START_REF] Rios | Derivative-free optimization: a review of algorithms and comparison of software implementations[END_REF], [Custódio, Scheinberg, and Nunes Vicente 2017]. The references for the respective solvers are: ASA [START_REF] Ingber | Very fast simulated re-annealing[END_REF]], BOBYQA [Powell 2009], CMA-ES [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF], COBYLA [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF]],

HOPSPACK [START_REF] Plantenga | Hopspack 2.0 user manual[END_REF]], IMFIL [START_REF] Kelley | Users Guide for IMFIL version 1.0[END_REF]], DAKOTA [START_REF] Adams | DAKOTA, a multilevel parallel object-oriented framework for design optimization, parameter estimation, uncertainty quantication, and sensitivity analysis: version 5.0 user's manual[END_REF], DFO [START_REF] Conn | A derivative free optimization algorithm in practice[END_REF], NEWUOA [START_REF] Powell | A direct search optimization method that models the objective and constraint functions by linear interpolation[END_REF]], NOMAD [START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear optimization with the MADS algorithm[END_REF],

PSwarm [START_REF] Vaz | PSwarm: a hybrid solver for linearly constrained global derivative-free optimization[END_REF], SID-PSM [START_REF] Custodio | SID-PSM: A pattern search method guided by simplex derivatives for use in derivative-free optimization[END_REF], SNOBFIT [START_REF] Huyer | SNOBFITstable noisy optimization by branch and t[END_REF] , TOMLAB [START_REF] Holmström | User's Guide for TOM-LAB 7[END_REF], GLODS [START_REF] Custódio | GLODS: global and local optimization using direct search[END_REF], PDF-MPC package [START_REF] Alamir | tiliztion of nonliner systems using reedingEhorizon ontrol shemesX prmetrized pproh for fst systems[END_REF]].

Fundamentals of Nonlinear Model Predictive Control

Introduction

The general formulation for a non-linear optimal control problem (OCP) is described by the following form min x(.),u(.)

T 0 L(x(t), u(t)) dt + ψ(x(T )) subject to ẋ = f (x(t), u(t)), ∀t ∈ [0, T ] g(x(t), u(t)) ≤ 0, ∀t ∈ [0, T ] h(x(t), u(t)) = 0, ∀t ∈ [0, T ] x(0) = x 0 , r(x(T )) ≤ 0 (2.5)
where x ∈ R n , u ∈ R m represents the state vector and input vector respectively, L : R n × R m → R and ψ : R n → R represents the stage cost and terminal cost respectively, ) is solved at every sampling period (τ ) and only the optimal input over one sampling period is injected into the system i.e. u * ([0, τ ]). In short, NMPC is non-linear OCP put in practice and it is also important to note that an implicit feedback is formed in the NMPC problem formulation as with receipt of new state information at every sampling period, the inputs are recomputed. Therefore, the optimal input is a function of the initial state as well (u * (x 0 )) [START_REF] Maciejowski | reditive ontrolX with onstrints[END_REF]].

f : R n × R m → R n , g : R n × R m → R p , h : R n × R m → R k
Fundamentally, there are three methods (illustrated in Fig. In this chapter, only the direct methods are described with considerable details. However, for the sake of completeness, a brief overview on the HJB equation/DP and PMP methods are given below.

HJB equation or DP method is based on the principle of optimality where the non-linear OCP is solved from the terminal (nal) cost to the initial cost in a backward fashion. The cost optimization method is divided into several stages (starting from terminal to initial stage costs) and the optimal input function or policy is computed for each stage recursively and nally, the optimal input policy is recovered for the initial stage and the initial state. It is important to note that HJB equation is used in the context of continuous time systems and DP method is used in the context of discrete time systems. The DP method in practice is typically solved by means of tabulation method and HJB equation is solved by solving a constrained Partial Dierential Equation (PDE). The major downside of the method is it suers from the curse of dimensionality and doesn't scale for higher dimensional systems [START_REF] Bertsekas | hynmi progrmming nd optiml ontrol[END_REF]. However, recently DP method has received considerable attention due to reinforcement and deep reinforcement learning community and the method in these elds are termed as approximate dynamic programming (ADP) or bootstrapping [START_REF] Bertsekas | einforement lerning nd optiml ontrol[END_REF]].

The Pontryagin minimum/maximum principle (PMP), also known as indirect method is based on the calculus of variations, where the objective is to minimize/maximize the Hamiltonian of the non-linear OCP. The method follows the philosophy of rst optimize and then discretize and is solved by means of gradient methods or shooting methods or collocation methods. The method can treat large scale systems and also solve boundary value problems (BVPs). Some of the major downsides include instability of the nonlinear ODE obtained for the state and co-state equations, requirement of explicit expression for the input and singular arcs are dicult to deal with and also, only necessary conditions for local optimality is obtained [Liberzon 2011]. However, the method is mostly used in the eld of aerospace engineering for computing the optimal orbit trajectory for satellites and spacecrafts. For more details, refer [START_REF] Naidu | yptiml ontrol systems[END_REF]].

Direct methods

The direct methods or direct transcription methods transcribes the problem into a nite dimensional problem, which yields an approximate solution to non-linear OCP (2.5) by virtue of NLP solvers. Direct methods follow the philosophy of rst discretize and then optimize.

The direct methods are classied into a) Single shooting methods, b) Multiple shooting methods and c) Direct collocation methods [START_REF] Rathai | GPU-Based Parameterized NMPC Scheme for Control of Half Car Vehicle With Semi-Active Suspension System[END_REF]. Single shooting

Optimal Control Problem

Hamilton-Jacobi-Bellman Equation Pontryagin Minimum/Maximum Principle Direct Methods

Single Shooting Methods

Multiple Shooting Methods

Collocation Methods The input u(t), ∀t ∈ [0, T ] is nitely parameterized by piecewise constant vector set U at an integer multiple of the sampling period (τ ) over the prediction horizon. With this representation, the input can be expressed with u(t) = µ(t; U), which is a piecewise continuous input signal.

The dynamics (dened by the ODE/DAE in (2.5)) of the system is numerically simulated for the given input signal u(t) = µ(t; U), and the solution is compactly expressed with x(t) = φ(t; U, x(0)), which is evaluated at N discrete time instants

T d = {t 0 , t 1 , . . . , t N } ⊂ [0, T ].
The time stamps T d typically corresponds to the integer multiple of the sampling period and it is utilized to discretize the dynamics, objective function and the constraint functions listed in (2.5). A ODE/DAE solver is utilized to simulate the system over the prediction horizon [START_REF] Ascher | gomputer methods for ordinry di'erentil equE tions nd di'erentilElgeri equtions[END_REF]]. With the aforementioned assumptions, the non-linear OCP is transcribed to a generic NLP framework, which is expressed with minimize w F (w)

subject to G(w) ≤ 0, H(w) = 0 (2.6)
where, w is the optimization variable which depends upon the direct method formulation, and F (w), G(w) and H(w) represents the objective, inequality constraints and equality constraints respectively. The transcribed NLP in equation ( 2.6) solves the non-linear OCP (2.5) and this is the crux of direct methods.

Remark -The ODE/DAE solver mentioned in the assumption is to be considered as a computer code and not in terms of algebraic equations. The ODE/DAE solver is a simulator which takes the numerical input trajectory u(t) = µ(t; U) and outputs the numerical state trajectory which is utilized in the objective and constraint functions for the NLP problem mentioned in (2.6). Thus, the objective and the constraint functions is embedded with the ODE/DAE solver code and is ought to be deemed as computer codes (function code). Under conditions of twice dierentiability of all the functions (codes) listed in (2.5), the Jacobians and Hessians for the NLP solver are numerically obtained by methods such as nite dierences, algorithmic dierentiation etc. [Griewank and Walther 2008] (also known as oracles in optimization parlance) and this information aids the optimization procedure.

Direct single shooting method

The direct single shooting method also known as sequential method eliminates the dynamics equality constraint in the non-linear OCP (2.5) by forward simulation and thus, removing the state variables from the OCP NMPC problem [START_REF] Rathai | GPU-Based Parameterized NMPC Scheme for Control of Half Car Vehicle With Semi-Active Suspension System[END_REF]. This reduces the optimization problem only to the input variables, which is obtained from the following NLP problem.

min

U N k=0 L(φ(t k ; U, x(0)), µ(t k ; U)) ∆t + ψ(φ(t N ; U, x(0))) subject to g(φ(t k ; U, x(0)), µ(t k ; U)) ≤ 0, ∀t k ∈ T d h(φ(t k ; U, x(0)), µ(t k ; U)) = 0, ∀t k ∈ T d r(φ(t N ; U, x(0))) ≤ 0 (2.7)
The objective is discretized by means of Riemann sum at time stamps T d and ∆t is the temporal dierence between t k+1 and t k , which usually coincides with the sampling period (τ ). The states are replaced with the ODE/DAE simulator evaluated at these time stamps in the objective as well as the constraints. The optimal solution U * obtained from (2.7) and the rst input U * (0) is injected into the system and the process is repeated in receding horizon manner. The prime benet of the method is, it only require few degrees of freedom, however

T 0 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 8 µ(t 0 ; U) µ(t 1 ; U) µ(t 2 ; U) µ(t 3 ; U) µ(t 4 ; U) µ(t 5 ; U) µ(t 6 ; U) µ(t 7 ; U) µ(t 8 ; U) µ(t 9 ; U) x(t) = φ(t; U, x(0))

Direct multiple shooting method

The direct multiple shooting method, also a part of simultaneous methods, retains the state variables as optimization variables and this increases the number of decision variables in the optimization formulation in the non-linear OCP (2.5). The ODE/DAE solver simulates the system over multiple time intervals i.e. [t k , t k+1 ], ∀k = {0, 1, . . . , N -1} simultaneously and the nal state value x(t k ) for the simulation in the each interval [t k , t k+1 ] is stipulated to obey the dynamics of the system, which is enforced by equality constraints [START_REF] Rathai | GPU-Based Parameterized NMPC Scheme for Control of Half Car Vehicle With Semi-Active Suspension System[END_REF]]. This method is implemented when the dynamics of the system (ODE/DAE equations) is numerically ill-conditioned and for larger prediction horizon. By dividing the prediction horizon into sub-intervals and simulating the system with respect to each subinterval, the stability property of the simulation is retained. The NLP optimization problem is formulated as min The optimal solution for the above optimization problem yields both the optimal state trajectory and optimal input sequence. As per the standard receding horizon policy, the rst input U * (0) is applied to the system and this procedure is repeated in the future. The However, the downside is that the optimization is carried out over an increased number of variables and a good initialization for the NLP solver is required for faster convergence to the optimal/suboptimal solution (this can be ameliorated by warm start procedure). A pictorial illustration of the method is shown in Fig. 2.4. For more details, refer [START_REF] Gros | From linear to nonlinear MPC: bridging the gap via the realtime iteration[END_REF]].

U ,{x(t 1 ),...x(t N )} N k=0 L(x(t k ), µ(t k ; U)) ∆t + ψ(x(t N )) subject to g(x(t k ), µ(t k ; U)) ≤ 0, ∀t k ∈ T d h(x(t k ), µ(t k ; U)) = 0, ∀t k ∈ T d x(t k+1 ) -φ(t k+1 ; x(t k ), µ(t k ; U)) = 0, ∀t k ∈ T d x(0) = x 0 , r(x(t N )) ≤ 0 (2.8) T 0 t 1 t 2 t 3 t 4 t 5 t 6 t 7 t 8 µ(t 0 ; U) µ(t 1 ; U) µ(t 2 ; U) µ(t 3 ; U) µ(t 4 ; U) µ(t 5 ; U) µ(t 6 ; U) µ(t 7 ; U) µ(t 8 ; U) µ(t 9 ; U) x(t k+1 ) = φ(t k ; U, x(t k ))

Direct collocation method

Direct collocation methods is a part of simultaneous methods, where the ODE/DAE simulator is expunged from the multiple shooting formulation (2.8) and is replaced with algebraic equal-Chapter 2. Theoretical background ity constraints enforced at the collocation points [START_REF] Rathai | GPU-Based Parameterized NMPC Scheme for Control of Half Car Vehicle With Semi-Active Suspension System[END_REF]. The solution for the dynamics (ODE/DAE) system is typically represented by means of orthogonal polynomials such as Lagrange polynomials and the innite ODE/DAE system is nitely discretized at the orthogonal collocation points dened by the polynomials. Therefore, the ODE/DAEs are gridded and transcribed into several equality constraints enforced at these collocation points and this method is also known as orthogonal collocation method. With this transformation, the NLP problem is expressed with min

U ,{x(t 1 ),...x(t N )} N k=0 L(x(t k ), µ(t k ; U)) ∆t + ψ(x(t N )) subject to g(x(t k ), µ(t k ; U)) ≤ 0, ∀t k ∈ T d h(x(t k ), µ(t k ; U)) = 0, ∀t k ∈ T d Ψ(x(t k+1 ), x(t k ), µ(t k ; U)) = 0, ∀t k ∈ T d x(0) = x 0 , r(x(t N )) ≤ 0 (2.9)
The fundamental dierence between (2.8) and (2.9) is Ψ, which denes the dynamics of the system at collocation points and also, unlike multiple shooting method, there are no online simulations of the ODE/DAE system. Direct collocation methods are typically suited for sti ODE/DAE systems. The benets of direct collocation methods are a) suitable for large scale problems and also, the NLP is sparse in nature, b) suitable for unstable systems, however the downside is that the method is not adaptive for new grid and changes in NLP dimensions.

For more details, refer [START_REF] Gros | From linear to nonlinear MPC: bridging the gap via the realtime iteration[END_REF]].

A Gentle

Introduction

As mentioned previously in Chapter 1, suspension systems play a vital role for control of vertical dynamics of vehicles for guaranteeing comfort and safety for the on-board passengers.

The seemingly simple task of control poses to be daunting under the presence of multiple nonlinearities, physical constraints and specication over objective requirements for the system.

Thereby, it is of paramount importance to account for these issues during control design for better eciency and prolonged endurance of the suspension system. Along with these control design requirements, it is also important to address the issue under the ambit of the EM-PHYSIS project [EMPHYSIS 2017], where the semi-active suspension system based vehicle models are exported from modeling tools such as AMESIM, OpenModelica, SimulationX etc.

as eFMU containers. Typically, the eFMU encapsulates the model of the system either as C source codes or binary executables and therefore, the generated eFMU could be considered as a black box model. Under such circumstances, where the knowledge of the system is completely obscured from the control engineer, the control design problem becomes more complicated in nature. Predicated upon these requirements, in this thesis a simulation-optimization based control strategy is proposed namely parameterized NMPC (pNMPC) scheme to handle this control problem. The pNMPC controller is tested, veried and validated for the INOVE test platform model by presuming the availability of the knowledge of the system, so that the method is scalable and copes up with black box model oriented eFMUs for future use case.

MPC is indisputably one of the most advanced and ecient control design methodology.

However, despite its enormous advantages in terms of optimal performance and constraint satisfaction, one of the major shortcoming is that the entire MPC controller hinges upon the model utilized in the control design. Due to the predictive nature of the controller, utilizing an erroneous model would ensue poor performance due to the mismatch of models between the plant and the controller. Thus, it is important to build a high delity model such that the mismatch is reduced and tangible performance benets from the MPC controller are obtained. Despite the mathematical model could be a black-box model or white-box model, the accuracy of the model with respect to the ground truth play a pivotal role in determining the performance of the MPC controller. Given this prelude and in the same spirit, in this chapter, the rst part addresses the problem of modeling and parameter identication for the ER semi-active suspension system for the INOVE test platform of a quarter car system. The obtained model could be conceptualized as a black-box model such as binary executables or eFMU/FMUs and the pNMPC scheme could be applied for this model. In the second part, the obtained model parameters are utilized for design and implementation of the proposed pNMPC scheme via HiL simulations on dSPACE MABXII and the INOVE test platform.

Related works

There have been several semi-active suspension system control design methods developed in the past such as Skyhook proposed in [START_REF] Karnopp | Vibration control using semiactive force generators[END_REF], Acceleration Driven Damping (ADD) proposed in [START_REF] Savaresi | Identication of semiphysical and black-box non-linear models: The case of MR-dampers for vehicles control[END_REF], Mixed Skyhook-ADD (SH-ADD) proposed in [START_REF] Savaresi | Mixed Sky-Hook and ADD: Approaching the Filtering Limits of a Semi-Active Suspension[END_REF], LPV/H-∞ based control methods proposed in [START_REF] Do | An LPV control approach for semiactive suspension control with actuator constraints[END_REF], [START_REF] Sename | Some LPV Approaches for Semi-active Suspension Control[END_REF]. A detailed literature review of dierent control strategies are presented in [START_REF] Tseng | State of the art survey: Active and semi-active suspension control[END_REF] and [START_REF] Poussot-Vassal | Survey on some automotive semi-active suspension control methods: A comparative study on a single-corner model[END_REF]]. Despite the several mentioned control strategies provide good performance, these methods adopt the state and input constraints in an ad-hoc fashion and not completely into control design. In [START_REF] Nguyen | LPV approaches for modelling and control of vehicle dynamics: application toa small car pilot plant with ER dampers[END_REF]], a robust control approach is applied by taking into account state and input constraints into control design, however, the method has several limitations and conservativeness. This exclusion of fully incorporating the constraints into control design might deteriorate the system performance and not fully utilize the potential of the semi-active damper system. To circumvent this problem, MPC based approach provides an elegant way of tackling the system constraints and objectives in control system design.

Another key advantage of MPC based approach is the ability to incorporate future road information (road preview) into control design which could improve the performance of the system in many folds.

Over the last decade, there has been several research contributions on MPC based approach for control of semi-active suspension systems. In [START_REF] Canale | Semi-active suspension control using "fast" model-predictive techniques[END_REF], a Fast MPC method is proposed where the optimal control input is computed oine by means of set membership approximation technique, however the model can not incorporate dynamic information into problem formulation such as road prole, variation in system parameters etc. In Hybrid MPC approach proposed in [START_REF] Giorgetti | Hybrid model predictive control application towards optimal semiactive suspension[END_REF]], the system is modeled as hybrid dynamical system and the optimal control input is computed oine by solving a multi-parametric program for a mixed-integer quadratic program (MIQP) and the method suers from similar shortcoming as for Fast MPC method. In [START_REF] Cseko | Analysis of the explicit model predictive control for semi-active suspension[END_REF], a detailed analysis of explicit MPC for semi-active suspension system is conducted.

In [START_REF] Gohrle | Active suspension controller using MPC based on a full-car model with preview information[END_REF]], a preview information based MPC scheme is proposed for control of suspension system for a full car model and the model is presumed to be a LTI system. In [Nguyen et al. 2016a], MPC for semi-active suspension system is implemented for full car model, where a MIQP problem is solved online, however, the sampling period is too high for practical implementation. In [START_REF] Morato | Design of a fast real-time LPV model predictive control system for semi-active suspension control of a full vehicle[END_REF]], a fast RT Linear Parameter Varying (LPV) MPC scheme is proposed for control of semi-active suspension system for a full vehicle to overcome the computational issues in [Nguyen et al. 2016a] by modeling the system by means of LPV model.

Chapter contributions

The main contribution of this chapter are System identication and parameter estimation of ER semi-active damper system for the INOVE test platform is explained in detail along with the design of experiments.

The obtained model parameters are utilized for implementation of the proposed pNMPC method for a quarter car system. It is also important to note that the estimated model and its parameters are utilized for all ER damper models throughout the thesis.

Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system

The working principle of pNMPC controller is explained in detail along with its application for control of semi-active suspension system. The design requirements such as objectives, constraints and the pNMPC OCP formulation for the suspension control system are laid down rmly.

The proposed method is implemented in HiL simulation on dSPACE MABX II and also, compared against a linearization based MPC using CVXGEN [Mattingley and

Boyd 2012] solver. The proposed method is experimentally implemented on the INOVE test platform and the performance was compared against other controllers.

Control oriented ER semi-active damper modeling and parameter identication

In general, any semi-active damper modeling exercise can be broadly classied into a) parametric and b) non-parametric based approach. Under the former regime, the structure of the model is dictated by the physics of the system, which is modeled by means of rst principles techniques and by contrast, the latter method obscures the underlying physics of the system and this leads to a exible model structure, which is modeled by means of empirical techniques.

[ [START_REF] Butz | Modelling and Simulation of Electro-and Magnetorheological Fluid Dampers[END_REF] provides a detailed survey on dierent types of parametric/nonparametric modeling for semi-active suspension systems. In this work, the parametric method is of primary interest and some of the popular parametric models (to name a few) include the Bingam model [START_REF] Stanway | Non-linear modelling of an electrorheological vibration damper[END_REF], the phenomenal Bouc-Wen model [START_REF] Spencer | Phenomenological model for magnetorheological dampers[END_REF]], the nonlinear viscoelastic-plastic model [START_REF] Kamath | A nonlinear viscoelastic-plastic model for electrorheological uids[END_REF],

the nonlinear bi-viscous model [START_REF] Stanway | Applications of electro-rheological uids in vibration control: a survey[END_REF]], Guo's damper model [START_REF] Guo | Dynamic modeling of magnetorheological damper behaviors[END_REF] etc. In this work, Guo's damper model is utilized to describe the ER semi-active suspension system for the INOVE test platform due to its simplicity and parsimonity to describe the damper characteristics.

Vehicle modeling -Quarter car model

The vertical dynamics model for the quarter car system equipped with ER semi-active damper system (around the equilibrium) is dened with the same quarter car model described in Chapter 1, Section 1.1 which is dened by

m s zs = -k s (z s -z us ) -u m us zus = k s (z s -z us ) + u -k t (z us -z r ) (3.1)
where, m s , m us are the sprung mass and unsprung mass respectively, k s , k us are the stiness coecients of suspension system and the tire respectively, z s , żs are the sprung mass position and velocity respectively, z us , żus are the unsprung mass position and velocity 3.3. Quasi-static nonlinear ER damper model 45 respectively, z r is the vertical road position or disturbance and u is the force exerted due to the ER semi-active damper system. The state vector is represented with x = [z s z us żs żus ] T .

Quasi-static nonlinear ER damper model

The ER semi-active damper force u is expressed using the quasi-static nonlinear damper model (Guo's damper model) [START_REF] Guo | Dynamic modeling of magnetorheological damper behaviors[END_REF] with

u = f c φtanh(a 1 żdef + a 2 z def ) + c 0 żdef (3.2)
where, c 0 , f c , a 1 and a 2 represents the viscous damping coecient, dynamic yield force of the uid, and hysteresis coecient of the damper model respectively. z def = z s -z us and żdef = żs -żus represents the deection position and velocity ∀i ∈ {l, r} between the chassis and wheel respectively. φ is the PWM-DC input for the system such that φ ∈ U, where U := [φ min , φ max ] such that 0 ≤ φ min < φ max ≤ 1. For reasons of safety and operational requirements, the minimum and maximum bounds for the PWM-DC signal were set to φ min = 0.1 and φ max = 0.35. To study the dynamic characteristics of the ER-SA damper system, it is important to estimate the response time of the system, which is indirectly estimated by nding the peak response time of the system (T r ). The following experiment was conducted to estimate T r :

The platform was excited with a pseudo binary random sequence (PRBS) road prole with an amplitude of 5 mm for a duration of 20 s, illustrated in Fig. 3.1.

The PWM-DC signal was ipped from φ min to φ max at time 10 s (a step change in input) and the ER-SA damper force (u) was measured, illustrated in Fig. 3.2.

In order to zero-in the point of transition (high frequency content), which provides the necessary cues for T r estimation, time-frequency analysis was performed by means of wavelet transform. The at-most period for PWM-DC transition to completely capture the dynamical behavior of the ER-SA damper system is ascertained.

The look ahead period or the prediction horizon for the MPC controller is computed.

Estimation of sampling time (T s ) for the damper system, which is computed using the general measure with T s ∈ [ T r 10 , T r 5 ] is computed [START_REF] Astrom | gomputer ontrolled systemsX theory nd design[END_REF]. This also accounts for all the delays in the system i.e. sensors and actuator delays. However, the natural sampling time of the system is 5 ms.

Design of experiments

In order to obtain the best model parameters for the ER semi-active damper system, it is imperative to conduct informative experiments and collect the input/output data that captures the dynamic behavior of the system. Conditioned upon the previous requirement, the test involved the following scenario:

PRBS signal based road excitation with an amplitude of 5 mm for a duration of 20 s.

A PRBS based input PWM-DC signal between the interval [φ min , φ max ] with a holding period of T r . semi-active suspension system

The rationale to adopt this scenario is to induce persistent excitation and minimize the crest factor for input design [START_REF] Ljung | ystem identi(tionX theory for the user[END_REF]]. The ER semi-active damper system was operated upon the aforementioned scenario and all the input/output data were collected for parameter estimation stage. For the purpose of illustration, Fig.

χ 2 (θ) = 1 N N i=1 D i Y -ψ(D i X , θ) 2 (3.3)
where, ψ is the estimated function for the data tting problem for the output dataset, which in this case is the ER semi-active suspension force, i.e. u. The optimal parameters are computed by solving the following nonlinear optimization problem θ * = argmin θ∈Θ χ 2 (θ) (3.4) where, Θ ⊂ R 4 is the constraint set for the parameters, which is dened with Θ

:= {θ ∈ R 4 | {θ 1 , θ 4 } ∈ R ≥0 , {θ 2 , θ 3 } ∈ R}.
The estimated model was validated using K-fold cross validation method with K = 5 and the accuracy of the model was estimated to 4.65 units.

The estimated model parameters are listed in the Table 3.1. Fig. 3.6 illustrates the predicted vs measured ER semi-active damper force (u) for a single track of input/output data. 50 Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system 3.5 pNMPC design requirements for semi-active suspension system

Objective requirements

The dichotomy of the objective design for the semi-active suspension system for a quarter car system could be both qualitatively and quantitatively classied into a) Comfort objective and b) Road Holding objective [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]].

Comfort objective (J com t ): Qualitatively, the prime goal of the comfort based objective design is to guarantee the comfort for the on-board passengers. The human body is sensitive to certain frequencies and it is of paramount importance to mitigate the eects of vibrations at these spots of the spectrum. Quantitatively, this tantamount to minimizing the vertical acceleration of the chassis (z s ). The comfort objective for the given look ahead period T l is expressed as

J com T l = T l 0 (z s (t)) 2 dt (3.5)
Road holding objective (J rh t ): Qualitatively, the prime goal of the road holding based objective design is to guarantee that the wheel is always in contact with the road. The requirement of this objective is crucial in control of longitudinal and lateral dynamics of the vehicle. Quantitatively, this objective corresponds to the requirement of minimizing the displacement between the road and the wheel (z us -z r ). The road holding objective for the given look ahead period T l is expressed as

J rh T l = T l 0 (z us (t) -z r (t)) 2 dt (3.6)
It is also important to note that both the objectives are conicting in nature. Thus, the objective for the semi-active suspension system holistically at the current time instant t is expressed as

J obj T l = θ 1 J com T l + θ 2 J rh T l (3.7)
3.6. Parameterized NMPC 51

Where, θ 1 and θ 2 are the weighting coecients between comfort and road holding objective and also, the coecients form a convex combination between the two objectives such that θ 1 + θ 2 = 1 and θ 1 , θ 2 ≥ 0.

Constraint requirements

The constraints for the semi-active suspension system primarily arises from the physical limitations of the system. These are hard constraints and must be handled systematically to prevent weariness of the system components. Dynamics constraint: The nonlinear equality constraints due to dynamics of the system dened in (3.1) and (3.2). It is important to note that this is eradicated via simulation aka direct single shooting method.

Road disturbance assumption: Under the consideration without road preview information, it is not uncommon to presume a constant road disturbance input measured at the current time instant t for the entire future horizon for the NMPC problem, i.e. d + = d(t).

Parameterized NMPC

The proposed parametrized NMPC approach is based on simulation methods, i.e. an explicit/implicit ODE solver is utilized to simulate the non-linear system in equation ( 3.1), (3.2) to determine the evolution of the states for a set of input sequences over the horizon. In the context of eFMI/FMU, the underlying system could be treated as black-box models and the solvers could include either Ordinary Dierential Equations (ODEs) or Dierential Algebraic Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system Equations (DAEs) solvers [START_REF] Ascher | gomputer methods for ordinry di'erentil equE tions nd di'erentilElgeri equtions[END_REF]], which could either be provided by the simulation environment or inbuilt within the eFMU/FMU containers. The optimal input sequence is elicited from the simulations which minimizes the objective function and satises the constraint requirements, which are handled algorithmically [Rathai et al. 2018]. The proposed parametrized NMPC algorithm is sequentially presented as follows Algorithm:

1. The input φ of the non-linear system in equation ( 3.1), (3.2) is nitely parameterized in time with N δ equidistant points over the look ahead period T l with {δ 0 . . . δ N δ -1 } time stamps with an interval of

T l N δ
and T l = δ N δ -1 and in space, the set U is discretized with N s points such that φ ∈ {φ 1 , . . . φ Ns } ⊂ U, where φ i is a discretization point in U. The input sequence over the horizon is compactly represented with µ(δ j |{φ i (δ j )} Ns i=1 , t), ∀j ∈ {0, . . . N δ -1}, i.e. at a given time instant δ j , there exists N s possible input values and this spans for all given time stamps.

2. The explicit/implicit ODE solver for the non-linear system in equation ( 3.1), (3.2) is simulated for all input sequences along space and time.

3. The optimal control sequence is computed with respect to the objective and constraints by plugging the simulated trajectory onto the cost function and the constraint functions.

The constraints are handled algorithmically that if a particular input sequence violates the constraints, then the input sequence is discarded and the solver is proceeded with another control sequence until the minimum cost is obtained. 4. In case, if no input sequence satises the constraints, then the input sequence which least violates the constraints is considered as the optimal input sequence.

5. This procedure is repeated in receding horizon policy method at every sampling period (T s ) and the optimal control input is φ * (0) = φ * (δ 0 ).

For the considered case of quarter car semi-active suspension system, N s is assumed as variable (space discretization) and N δ = 1 (time discretization) and the solver utilized is a simple fourth order explicit Runge-Kutta (RK) method with xed integration step h = 1 ms.

The model parameters for INOVE quarter car platform and proposed MPC design are listed in Table 3.2. It is important to note that the ascertained sampling period from the system identication tests is T sest = 0.023 and the natural sampling period of the system or the sampling period of the DAQ is T s = 0.005. For experiments on the test bench the former was used and for tests on dSPACE MABXII the later was used. The proposed pNMPC controller was compared against a linearization based MPC controller.

The performance and computation time were compared with each other to gauge the potential of the proposed pNMPC method.

Linearization based MPC design

The fundamental assumption for the linearization based MPC design is to linearize the nonlinearities present in the system (i.e. constraints and dynamics) by means of rst order Taylor series expansion and then, the problem is casted as a linear MPC problem, i.e. a convex QP problem. This procedure is repeated at every operating point and a linear MPC is solved at every operating point in receding horizon fashion. The rst order linearization of the quasistatic nonlinear damper model (3.2) under the modied input at a given operating point

P i = (x i , φ i ) is expressed as u P i (x k , φ k ) = u(P i ) + ∇ φ u| (P i ) ∆φ + ∇ x u| (P i ) ∆x (3.8)
Where, ∆x = x k -x i and ∆φ = φ k -φ i are the state and input deviation variables with respect to the operating point P i . The rst order linearization of the nonlinear dynamics at the operating point P i in continuous time is expressed as

∆ ẋ(t) = A i c ∆x(t) + B i c ∆φ(t) + B cd d(t) (3.9)
Where, A i c ∈ R 4×4 and B i c ∈ R 4×1 are the linearized system and input matrices at P i . The obtained continuous time matrices are converted to discrete time matrices by means of zero 54 Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system order hold (ZOH) method with a sample time T s . The discrete-time linearized state space equation at the point P i is expressed as

∆x + = A i d ∆x(k) + B i d ∆φ(k) + B dd d(k) (3.10) Where, A i d ∈ R 4×4 , B i d ∈ R 4×1 and B dd ∈ R 4×1
are the discrete-time system matrix, input matrix and disturbance matrix. The linearized MPC nite time optimal control problem (FTOCP) at the point P i with with x 0 = x(0) and with d 0 = d(0) is casted as a convex QP problem with horizon length N corresponding to the look ahead period T l which is described as

J * P i = min φ 0:N -1 ,x 1:N N -1 k=0 J obj k s.t. z s -z us ∈ [z min def , z max def ], ∀k = 1 . . . N d k+1 = d k , ∀k = 0 . . . N -1 u P i (x k , φ k ) ∈ [u, u], ∀k = 0 . . . N -1 φ k ∈ U, ∀k = 0 . . . N -1 (3.10), ∀k = 0 . . . N -1 (3.11)
Where, J * P i is the optimal objective function. The optimal control input at point P i to the actual system is φ * (0) = φ * 0 and this procedure is repeated in receding horizon policy method.

For initialization of the input for the next linearization point, the solution of the previous program of equation (3.11) is utilized, i.e. P i+1 = (x i+1 , φ * 1 ), where φ * 1 is the solution at time step 1 at P i point. The linearization is performed by precomputing the Jacobians a priori and is evaluated at every time instant.

Simulation analysis for pNMPC method

A detailed analysis is conducted for the proposed parameterized NMPC method for the quarter car semi-active suspension system for dierent cases. The parameterized NMPC method is simulated in MATLAB/Simulink environment and its closed-loop performance characteristics are investigated for dierent complexity factors i.e. dierent space discretization points (N s ) and computational time i.e. the control update period (τ φ ). The considered acid test is for the following scenario A chirp road prole with a frequency sweep between 0.1 Hz to 25 Hz with an amplitude of 1 mm for a duration of 10 s.

The control objective is selected to be comfort oriented design i.e. θ 1 = 1 (θ 2 = 0).

The control update period is a variable which is expressed with τ φ = γN s T s , where γ is the computational scale factor and T s = 0.5 × 10 -4 s.

The rationale behind this heuristic and the analysis is to comprehend the behavior of the proposed controller when executed in dierent computational resources for dierent complexities (N s ), control update period (τ φ ) and also, the analysis provides insight over the optimal s to be utilized for a given computational resource. Fig. 3.7 illustrates the normalized closed loop performance of the system for complexity factor (N s ) vs normalized closed loop objective (J norm CL ) for dierent computational scale factor (γ). The normalized closed loop objective (J norm CL ) is dened with respect to the objective of nominal passive sus- pension system dened as

J norm CL = J obj CL J pass CL (3.12)
where, J pass CL corresponds to the the closed loop objective for the nominal passive suspension system and J obj CL corresponds to the closed loop objective of the parameterized NMPC method. The curves in plot Fig. 3.7 illustrates the fact that the normalized closed loop objective (J norm CL )

for a given computational scale factor (γ) declines as the complexity factor gradually increases, however as the complexity factor increases more than a certain threshold, the normalized closed Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system loop objective (J norm CL ) increases due to the fact that the computational load is elevated and consequently, the control update period (τ φ ) is increased, which results in poor performance of the controller. The abscissa of the optimal point for the curves indicates the best/optimal complexity factor N * s for a given computational resource or computational scale factor (γ).

Remark -The proposed parametrized NMPC is of high interest for practical applications for a large set of semi-active dampers. Indeed it is worth noting that N s denes the set of damping coecients than can be used in real-time control. When that N s tends to innity this corresponds to a continuously-variable damper. When N s = 2 this corresponds to a 2states damper or to a min-max suspension control approach (as for SkyHook, and ADD and SH-ADD methods). This oers a large exibility for the implementation of several control methods for dierent damper types.

Computational eciency test

To test the computational eciency of the two methods, the maximum execution time is recorded for dierent complexity parameter of the controller. The complexity parameter for the linearization based MPC is selected to be the number of Newton steps/iterations for the QP solver and for the proposed approach, the number of space discretization points (N s ) is considered. The road prole is a chirp signal with amplitude of 1 mm and frequency sweep from 5 Hz to 25 Hz with comfort objective i.e. θ 1 = 1 (also applies for road holding objective θ 2 = 0). The test involves a chirp road prole (shown in Fig. 3.11) with amplitude of 1 mm and frequency sweep from 5 Hz to 25 Hz with comfort objective i.e. θ 1 = 1 (also applies for road holding objective i.e. θ 2 = 0). The complexity factor for linearization based MPC, i.e. number of iterations is 11 and for parameterized NMPC N s is 20. The root mean square (RMS) values for the simulations are listed in Table 3.3. For the comfort objective, the chassis acceleration is shown in Fig. 3.12. The dissipativity constraint due to non-linear modeling of ER damper is shown in Fig. 3.13. The results demonstrate better performance of the proposed approach compared to linearization based MPC for the considered scenario in RT considerations.

Comparison controllers

Modied Skyhook controller -Skyhook controller is one of the most prominent and well known controller for semi-active suspension system [START_REF] Karnopp | Vibration control using semiactive force generators[END_REF]. The modied skyhook controller is an extension to the skyhook controller where the ER-SA damper system's PWM-DC signal swings between minimum and maximum value conditioned upon a switch condition. Mathematically, the controller is expressed with

φ = φ max , if żs żdef ≥ 0 φ min , if żs żdef < 0 (3.13)
Nominal passive suspension -The nominal passive suspension is a typical passive suspension system, however the term nominal indicates that the PWM-DC for the ER-SA damper system is xed to the mean value of the minimum and maximum values of the PWM-DC, i.e. φ nom = φ min + φ max 2

. The value is held constant over the entire period of operation, which is 0.225.

Results and Implementation

The

Conclusions

This chapter has presented the pNMPC scheme for control of vertical dynamics of vehicle via ER semi-active suspension system for a quarter car model. The method was tested through HiL simulations as well as the INOVE test platform to validate and verify the performance with respect to the objective and constraints requirements. The method was compared against Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system the dynamical models are exported as eFMU containers, where at times, the knowledge of the system is not explicitly provided to the end user. The proposed pNMPC method does not require any internal details of the system rather only the model ought to be simulated to obtain the stae trajectory of the system. Thus, the proposed pNMPC method is scalable with eFMU containers and also, can cope up with black-box models. 66 Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active suspension system

Introduction

This chapter introduces the parallelized pNMPC scheme for control of semi-active suspension system for a half car vehicle. The method taps the power of GPU computing for parallelizing the pNMPC control scheme across several GPU multi-core processors. The method is generic in nature, however in thesis the proposed method is utilized for solving the automotive suspension control problem. The chapter is divided into two parts, where in the rst part the road model is not accounted into the OCP formulation and in the second part the road prole is modeled via a stochastic process and the resulting method is termed scenario stochastic pNMPC scheme (SS-pNMPC).

Some of the key aspects of the proposed methods are

Automotive standard -The method is amenable with the existing automotive standards as only simple math operations are required and also, independent of any back-end optimization solver.

Modeling exibility -The proposed SS-pNMPC method requires no a-priori assumptions on the system such as linear time invariant (LTI) dynamics, stationary property, Gaussian distribution of noise etc. This feature is highly sought these days as proposed in [START_REF] Guanetti | Stochastic MPC for cloud-aided suspension control[END_REF], the road model stochastic information can be dynamically obtained from cloud servers and this could be easily embedded in the OCP formulation.

Black-box model compatibility -When dealing with black box models such as FMUs [START_REF] Blockwitz | Functional mockup interface 2.0: The standard for tool independent exchange of simulation models[END_REF] or eFMUs, the model implementation details such as dynamical equations, system parameters etc. are concealed away from the end user due to protection of intellectual property (IP) rights. In the such a situation, designing a model based controller poses to be a serious challenge, however the proposed pNMPC controller is based on a simulation-optimization method and this obviates any need to stipulate the structure of dynamics, objective or constraint functions. The proposed approach is implemented on GPU to increase the computation throughput for the controller.

Ecient RT operability -The fast computation of the control input by means of GPU renders the method to be RT operable, especially for fast systems.

Validation/Verication on Embedded platforms -The proposed SS-pNMPC method was tested on multiple GPU based embedded boards to calculate the computation time and also to assess the RT feasibilty of the method. The method was tested on the NVIDIA Jetson embedded boards -Nano, TX1, TX2 and Xavier and the results looks promising for RT implementation.

Related works

There has been several research conducted in a recent past on solving the MPC problem which harnesses the potential of GPUs. A brief tutorial on dierent parallel architectures for MPC is proposed in [START_REF] Koegel | Parallel architectures for model predictive control[END_REF]. In [Gade-Nielsen, Dammann, and Jørgensen 2014],

several GPU based interior point methods were developed for linear MPC framework. There have been several research contributions on utilization of GPUs for solving stochastic control problems and in [START_REF] Abughalieh | A Survey of Parallel Implementations for Model Predictive Control[END_REF], a detailed survey on various types of parallel implementation of MPC methods are described. In [START_REF] Sampathirao | GPU-accelerated stochastic predictive control of drinking water networks[END_REF]], a scenario based Stochastic MPC (SMPC) method is proposed where the structure of the system is exploited and the problem is solved using proximal gradient method which is parallelized on GPU. In [Ohyama and Date 2017], a sampling based parallelized nonlinear MPC (NMPC) scheme is proposed and experimentally validated for control of inverted pendulum system. In [START_REF] Williams | Aggressive driving with model predictive path integral control[END_REF]], a path integral based MPC method is proposed and experimentally validated for a Remote Controlled (RC) car, the paper derives an input update rule for a stochastic optimal control problem based on information theoretic concepts and at every sampling period, multiple random scenarios are generated to update the input sequence. In [START_REF] Rogers | GPU-enabled projectile guidance for impact area constraints[END_REF]], a guidance law for guided projectiles is proposed, where the GPUs are utilized to generate RT scenarios to predict the impact point and probability of violating impact area constraints. In [START_REF] Hyatt | Real-time nonlinear model predictive control of robots using a graphics processing unit[END_REF], a GPU based RT NMPC scheme for control of robots is proposed where the method is termed as Nonlinear Evolutionary MPC (NEMPC) and is practically implemented for control of a 24 state pneumatically actuated continuum soft robot.

Concerning the control of semi-active suspensions system, as mentioned previously in Chapter 3, [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]] provides a comprehensive collection of all classical and modern control methods such as (to name a few) Skyhook, SH-ADD, Hybrid MPC with preview, H ∞ and LPV methods etc. However, in the line of research of application of Stochastic MPC for control of suspension systems, not many research have been conducted in the past. To the best of knowledge of the authors, in [START_REF] Guanetti | Stochastic MPC for cloud-aided suspension control[END_REF]] a cloud aided SMPC method is proposed for control of active suspension system for a quarter car vehicle. However, the method is proposed only for LTI systems and also, the second order conic program (SOCP) solver utilized in the method is computationally not tractable within the prescribed sampling period for higher order systems such as half/full car.

Chapter contribution

The main contributions of this chapter are

The pNMPC method proposed in Chapter 3 is augmented by parallelizing the method over the several multi-core streaming processors in the GPUs. By parallelizing the pN-MPC method over several multi-core processors, one can attain signicant reduction in computation time due to several simultaneous simulations of the dynamical system under several control congurations. The proposed parallelized pNMPC method is applied for control of semi-active suspension system for a half-car model with two congurations.

In the rst conguration, the road model is not included into the control design and the proposed controller is simply termed as parallelized pNMPC scheme. In the second conguration, the road model is modeled via a stochastic model using the ISO road pro-68 Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active suspension system le standand and the proposed controller is termed Scenario-Stochastic pNMPC scheme (SS-pNMPC).

The proposed methods were tested in simulation as well as in NVIDIA embedded boards to test, verify and validate the performance both in simulation as well as in RT scenario.

The methods were developed using C++ with CUDA-C compiler support and MAT-LAB/Simulink environment.

           m s zs = -i∈{l,r} F s,i I x θ = (l l F s,l -l r F s,r )
m us,l zus,l = (-F s,l + F t,l ) m us,r zus,r = (-F s,r + F t,r ) The equations of the model is expressed as follows

                           m s zs = -(F s,l + F s,r ) I x θ = (l l F s,l -l r F s,r ) m us,l zus,l = (-F s,l + F t,l )
m us,r zus,r = (-F s,r + F t,r ) żr,l = -αv x z r,l + ξ l żr,r = -αv x z r,r + ξ r vx = a x

H(jω)Ψ w H * (jω) = 2αv x σ 2 (αv x + jω)(αv x -jω) (4.3)
where, ω is the excitation frequency of the road input (in rad/s), Ψ w = 2αv x σ 2 represents the spectral density of the Gaussian white noise and σ 2 denotes the road roughness variance. H(jω) = 1 (αvx+jω) is the rst order shaping lter. This frequency response function can be recasted as a time-varying rst order auto regressive process, which is expressed with żr (t) = -αv x (t)z r (t) + ξ(t) (4.4) where, z r (t) is the road prole as dened in (4.2) and road disturbance drawn from a normal distribution dened with ξ(t) ∼ N (0, Ψ zr (t)). The parameters for dierent road surfaces are listed in Table 4.1. It is also important to note that the road roughness variance listed in Table 4.1 is scaled accordingly to suit the INOVE test platform.

Mathematical terminology

z s , θ represents the heave/chassis position and roll angle of the vehicle w.r.t. the centre of gravity (COG) respectively. z us,i , z r,i , ξ i ∀i ∈ {l, r} represents the wheel/unsprung mass position, vertical road prole of the vehicle and random disturbances respectively [START_REF] Guanetti | Stochastic MPC for cloud-aided suspension control[END_REF]. v x and a x denotes the longitudinal velocity and acceleration of the vehicle. m s , m us,l , m us,r represents the chassis mass, unsprung masses for the left and right corners. I x represents the moment of inertia along the roll axis. l l and l r represents the length of the chassis from the left and right corners with respect to COG. α represents the ISO road prole parameter. F s,i represents the chassis forces and F t,i represents the wheel forces ∀i ∈ {l, r} which are expressed with 5) suspension system where, k s,i and k t,i represents the stiness coecent of the SA suspension system and wheel respectively. z r,i and z us,i represents the vertical road displacement and unsprung mass position ∀i ∈ {l, r}. u i represents the actuation force obtained from the nonlinear SA damper model (see Section 4.2.5). z s,i represents the sprung mass displacement at each corner which are obtained from the following equations z s,l = z s + l l sinθ z s,r = z s -l r sinθ (4.6)

F s,i = -k s,i (z s,i -z us,i ) + u i F t,i = -k t,i (z us,i -z r,i ) (4.

Nonlinear quasi-static SA damper model

The SA damper force u i (4.4) is expressed by utilizing the Guo's damper force model [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF] ∀i ∈ {l, r} with

u i = c 0 żd,i + f c φ i tanh(a 1 żd,i + a 2 z d,i ) (4.7)
where k 0 , c 0 , f c , a 1 and a 2 represent the damper stiness coecient, viscous damping coecient, dynamic yield force of the uid, hysteresis coecient due to velocity and position respectively. φ i , ∀i ∈ {l, r} represents the duty cycle (PWM-DC) input signal which manipulates the damper characteristics online by changing the input voltage. z d,i = z s,i -z us,i and żd,i = żs,i -żus,i represents the deection position and velocity ∀i ∈ {l, r} between the chassis and wheel respectively. The values for the parameters are listed in Chapter 3, Table 3

= [z r,l , z r,r ] denote the disturbance vector. X ∈ R 8 , U ∈ R 2 and D ∈ R 2 ,
then the half car model in equation ( 4.1) can be compactly expressed with 3. Road disturbance assumption: The road disturbance is assumed to be constant over the prediction horizon (T l ) with D(0) = D 0 i.e. the road prole measured at the current time instant by means of observers [START_REF] Doumiati | Road prole estimation using an adaptive YoulaKu£era parametric observer: Comparison to real prolers[END_REF]]. There is no loss of generality with this formulation as it can be easily extendable with any road models such as ISO road proles or from road preview sensors.

Ẋ(t) = f (X(t), U(t), D(t))
The mixed input-state constraint set is compactly expressed with (X, U) ∈ Ω (X,U) ⊂ R 8 × R 2 , the input constraint set is compactly expressed with U ∈ Ω U ⊂ R 2 and the state constraint set is compactly expressed with X ∈ Ω X ⊂ R 8 .

MPC problem formulation

In summary, with the proposed objective and constraints function, the NMPC OCP is casted as J obj (X 0 , Γ 0 , D 0 , U(.)) = min U(.)

Γ 1 0 J acc T l + Γ 2 0 J roll T l subject to Ẋ(t) = f (X(t), U(t), D(t)) X 0 = X(0), D(t) = D 0 X(t) ∈ Ω X , U(t) ∈ Ω U (X(t), U(t)) ∈ Ω (X,U) (4.11)
where Γ 0 = [Γ 1 0 , Γ 2 0 ] with Γ 1 0 and Γ 2 0 being the convex combination weights between the two objectives i.e. Γ 1 0 + Γ 2 0 = 1. Once the optimal input trajectory is computed, the rst control action in U * is injected into the system and this procedure is repeated in receding horizon fashion. In this work, a constant input prole is assumed over the control horizon.

Parallelized pNMPC Method

The crux of the method is to nitely parameterize the input constraint set Ω U and by virtue of parallel computing methods, the system in equation ( 4.8) is simulated for every parameterized input given the current state, disturbance and objective specication information. The optimal input is elicited which minimizes the objective and satises the constraint in equation (4.11) (see [Rathai et al. 2018] for more details). The pseudo-code for the implementation of parallelized pNMPC is shown in Implementation 1.

The pseudo-code deserves some explanation to elucidate its working principle. The details enclosed in this section provides the reader with the need to know basics to handle parallel computing using the NVIDIA CUDA GPUs for implementation of the proposed parallelized pNMPC method. The pseudo-code is written in a manner to benet the reader to easily connect with CUDA-C implementation. For the respective application programming interfaces (APIs) and other syntax details refer to [START_REF] Sanders | ghe y exmpleX n introdution to generlE purpose q progrmming[END_REF].

Explanation of Implementation:

1. Initialization, I/Os and Syntax declaration: (b) The input variables for the method are X 0 , D 0 , Γ 0 and the output variable is U * . (c) The line 13 is a check condition to make sure the thread access is not exceeded.

(d) The lines 16 and 17 obtains the input combination and objective value from the functions GRID2D and ODESIM2OBJ. KERNEL {gs,bs} (d obj , d inp , X 0 , D 0 , Γ 0 )

d inp [r] ← GRID2D(i, j) 17: d obj [r] ← ODESIM2OBJ(X 0 , D 0 , Γ 0 , i,
φ i,l = φ min,l + i n φ l -1 (φ max,l -φ min,l ) 21:
φ j,r = φ min,r + j n φr -1 (φ max,r -φ min,r )

22:

return {φ i,l , φ j,r } 23: end function 24: function __device__ ODESIM2OBJ(X 0 , D 0 , Γ 0 , i, j) 25:

Obj = 0, Con = 0 26: U i,j ← GRID2D(i, j) 27:
for t loop = 0 : h : T l do 28:

X + ← X 0 + hf (X 0 , U i,j , D 0 ) 29: if X + / ∈ Ω (X,U) ∨ X + / ∈ Ω (X) then 30: Obj = MAX; Con = Con + N(Ω (X,U) , Ω (X) ) 31:
else 32:

Obj = Obj + h(Γ 1 0 J acc l l + Γ 2 0 J roll l l ) 33: end if 34: X 0 = X + 35:
end for 36:

return {Obj + Con} 37: end function

Analysis and Simulation results

The conducted simulation study can be broadly classied into three parts which are 1. Computational time (CT) analysis b/w CPU (serial) and GPU (parallel) for the proposed pNMPC method.

2. Comparative analysis b/w ACADO-qpOASES NMPC controller [START_REF] Houska | ACADO toolkit user's manual[END_REF] and the proposed parallelized pNMPC method.

3. Performance analysis with a road prole test. Both methods were implemented in MAT-LAB/Simulink on a Intel Core i7 PC and NVIDIA GTX 1050Ti with 768 CUDA cores Code generation option was utilized for ACADO-qpOASES NMPC controller and the proposed parallelized pNMPC was programmed in CUDA C and patched into Simulink with S-function.

Computational time analysis b/w CPU and GPU

The raison d'être for conducting this analysis is to emphasize the signicance of GPUs for solving huge simulations and viability of the approach for control of real SA suspension system. From Fig. 4.1, it is evident that the pNMPC method fares well in GPU compared to CPU in terms of mean CT per sampling period. The abscissa indicates the number of discretized inputs for the PWM-DC signal for both left (n φ l ) and right (n φr ) corner of the vehicle (i.e. n φ l × n φr input combinations). The road prole involved a chirp signal (same prole for both corners) with an amplitude of 2.5 mm and frequency sweep from 1 -14 Hz for a duration of 10 s. The objective was comfort (Γ 1 0 = 1, Γ 2 0 = 0).

Comparative analysis

A comparative analysis was conducted b/w ACADO-qpOASES NMPC controller and the proposed parallelized pNMPC method. The basis for conducting this study was to analyze computational time (CT), normalized closed-loop objective (NCLO), i.e. with respect to nominal damping (see Section 4.4.3) and feasibility analysis (FA), i.e. constraint satisfaction.

The study involved a gradual increase in the complexity parameter and the aforementioned criteria were recorded. The complexity parameter for ACADO-qpOASES NMPC controller was the number of Newton iterations (N s ) i.e. IMPLICIT_INTEGRATOR_NUM_ITS (see [START_REF] Houska | ACADO toolkit user's manual[END_REF]) and other settings were set to default and the number of discretization points {n φ l , n φr } (left/right corner of the vehicle) for parallelized pNMPC method. The road prole utilized was the same as mentioned in Section 4.4.1 and the objective was comfort (Γ

1 0 = 1, Γ 2 0 = 0).
The and indicates the infeasibility and feasibilty of the imposed constraints for the system. The recorded readings are listed in Discretized inputs (Complexity)

Road profile for ride handling

Left road profile Right road profile Roll angle (rad)

Ride handling performance

Parallelized pNMPC method ACADO-qpOASES Nominal damping in Fig. 4.6, assumes the inputs to be in continuum, which in practice might not render the required performance for the system. On the contrary, this can be easily included into the proposed parallelized parameterized NMPC (pNMPC) method and this is illustrated in Fig. 4.5.

4.5 Scenario-stochastic pNMPC scheme for control of semiactive suspension system 4.5.1 Mathematical model notations Let X = [z s , θ, z us,l , z us,r , żs , θ, żus,l , żus,r , z r,l , z r,r , v x ] denotes the state vector, U = [φ l , φ r ] denotes the input vector and Ξ = [ξ l , ξ r ] denote the disturbance vector. X ∈ R 11 , U ∈ R 2 and Ξ ∈ R 2 , then the half car model in equation ( 4.2) can be compactly expressed with

Ẋ(t) = f (X(t), U(t), Ξ(t), a x (t)) (4.
12)

It is important to note that all the state variables are assumed to be measured and this also includes the road prole [z r,l , z r,r ] at every sampling period.

SS-pNMPC design requirements 4.5.2.1 Objective requirements for SS-pNMPC controller

The objective design in time domain can be briey classied as a) comfort and b) ride handling objective [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]].

1. Comfort objective: The goal of the comfort based objective is to minimize the vertical acceleration of the chassis (z s ), governed by equation (4.2). The comfort objective for a given look ahead period T l is expressed with J acc T l (X(.), U(.), Ξ(.), a x (.)) =

T l 0 (z s (t)) 2 dt (4.13) 2.
Ride handling objective: The goal of the ride handling objective is to minimize the roll angle (θ) of the vehicle. The ride handling objective for a given look ahead period

T l is expressed with J roll T l (X(.), U(.), Ξ(.), a x (.)) = T l 0 (θ(t)) 2 dt (4.14)

Constraint requirements for SS-pNMPC controller

The constraints for the SA suspension system primarily arise from the physical limitations and secondarily from performance requirements [START_REF] Baumal | Application of genetic algorithms to the design optimization of an active vehicle suspension system[END_REF]]. For the MPC design, the included constraints are (b) Wheel rebound constraint (Performance): This bounds the deection position between the wheel and road. This is to ensure the tyre deection forces are bounded i.e. z us,i -z r,i ∈ [z rebmin,i , z rebmax,i ], ∀i ∈ {l, r}.

3. Road variable assumption: The vertical road displacement at the current time instant is assumed to be measured by means of adaptive road prole observers [START_REF] Doumiati | Road prole estimation using an adaptive YoulaKu£era parametric observer: Comparison to real prolers[END_REF] or from cloud servers [START_REF] Zhang | Cloud-Aided State Estimation of A Full-Car Semi-Active Suspension System[END_REF]].

Longitudinal acceleration assumption:

The longitudinal acceleration a x is assumed to be constant over the prediction horizon (T l ). In a real vehicle setting a x is typically obtained from Inertial Measurement Unit (IMU) of the vehicle.

The mixed input-state constraint set is compactly expressed with (X, U) ∈ Ω (X,U) ⊂ R 11 × R 2 and the input constraint set is compactly expressed with U ∈ Ω U ⊂ R 2 . The disturbance has a probabilistic support P which is normally distributed with Ξ ∼ N (0, Σ(X(t))), where Σ(X(t)) = diag(Ψ z r,l (X(t)), Ψ zr,r (X(t))). The variance is dependent on longitudinal velocity (v x ) of the vehicle, which is a state variable of the system (see Section 4.2.3).

SNMPC problem formulation

In summary, with the proposed objective and constraints function, the SNMPC OCP is casted as

J * obj (X 0 , Γ 0 , a x,0 ,U * (.)) = min U(.) E(Γ 1 0 J acc T l + Γ 2 0 J roll T l ) s.t. Ẋ(t) = f (X(t), U(t), Ξ(t), a x (t)) X(0) = X 0 , a x (.) = a x,0 , U(.) ∈ Ω U P[(X(.), U(.)) / ∈ Ω (X,U) ] ≤ η (4.15)
where X 0 represents the initial state vector and a x,0 represents the constant longitudinal acceleration over the prediction horizon. Γ 0 = [Γ 1 0 , Γ 2 0 ] with Γ 1 0 and Γ 2 0 being the convex combination weights between the two objectives comfort and ride handling respectively. The stochastic measure adopted for the total objective is the expectation operator (E). The mixed state-input constraints (Ω (X,U) ) are encapsulated in a probabilistic framework with a nite level of violation η 1. Once the optimal input trajectory is computed, the rst control action is injected into the system and this procedure is repeated in receding horizon fashion.

Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active suspension system

It is also important to note that the two objectives are conicting in nature (See [START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]). In this work, a constant control input prole is assumed over the control horizon.

As the semi-active suspension system is inherently stable the foregoing assumption is apposite for performance requirements for fast sampled systems. .

Method description

The method is an extension of the work proposed in [Rathai et al. 2018] and [Rathai, Sename, and Alamir 2019], where the stochastic road model is accounted in the model dynamics. The fundamental idea of the method is to parameterize the input set (Ω U ) into nite number of control inputs (similar to nite control set MPC) and for each parameterized control input, the SMPC is solved by performing MC simulations with several scenarios of the road prole.

From the simulations, the expected objective function is numerically obtained by empirical mean and a probabilistic constraint violation certicate (PCVC) is numerically obtained by computing the ratio between constraint violation and scenarios generated. The optimal input is selected by nding the minimum expected objective along with the consideration of PCVC less than or equal to the specied level (η). If none of the input satises the above criteria, then the input with the least PCVC is selected. The pseudo-code for the method is shown in Algorithm 2. The algorithm is explained in the following part. It is important to note is that the pseudo-code is an abstraction of the parallel programming paradigm. Several threads are spawned in parallel to achieve this task (See [START_REF] Sanders | ghe y exmpleX n introdution to generlE purpose q progrmming[END_REF]).

Explanation of Implementation:

1. Initialization, I/Os and Syntax declaration:

(a) The data initialization step sets the parameter values for the half-car model and the constraints.

(b) The input variables are X 0 , Γ 0 , a x,0 , n g , γ, where n g and γ are the number of input parameterization and number of scenarios respectively. The output variable is U i * , i.e. the optimal input vector injected into the system.

(c) The qualiers __CPU__ and __GPU__ denotes the function operation in CPU and GPU respectively. The entry point is SSpNMPC function.

SSpNMPC function:

(a) In line 2, the input set Ω U is nitely discretized into n g points and collected in the grid U 1:ng .

(b) From line 3-8, the parfor i.e. parallel for function is utilized to dispatch each and every discretized input from U 1:ng to GPU and for each input the SIM function is utilized to conduct N s number of MC simulations. The respective objective function Obj i [l] and constraint violation CV i [l] for the l th MC simulation and i th input are 4.5. Scenario-stochastic pNMPC scheme for control of semi-active suspension system 83 Implementation 2 SS-pNMPC pseudocode Data Initialization: Model/Constraint parameters Input: X 0 , Γ 0 , a x,0 , n g , γ

:= N s × N τ Output: U i * 1: function __CPU__ SSpNMPC(X 0 , Γ 0 , a x,0 ) 2:
U 1:ng ← grid(Ω U , n g ) 3:

parfor i ← 1 : n g do 4:

parfor l ← 1 : N s do 5: return U i * 16: end function 17: function __GPU__ SIM(X 0 , Γ 0 , a x,0 , U i ) 18: (c) The lines 9-15, looks for the optimal input i * index with minimum expected objective i.e. EObj i * and also satises the violation condition PCVC i * ≤ η. If the aforementioned condition is true, i * is obtained from the function indexmin over the vector EObj, otherwise, i * with the least violating constraint is obtained from the function indexmin over the vector PCVC. Once the index for the optimal input (i * ) is obtained, U i * is injected to the system. It is important to note that the optimal objective (J * obj ) and input (U * (.)) in (4.15) are EObj i * and U i * .

(Obj i [l], CV i [l]) ← SIM(X 0 , Γ 0 , a x,0 , U i ) 6: end parfor 7: EObj i ← l Obj i [l] Ns ; PCVC i ← l CV i [l] Ns 8: end parfor 9: if (EObj i * ≤ EObj ∀i\{i * } & PCVC i * ≤ η) then 10: i * ← indexmin(EObj) 11: J * obj ← EObj i * 12:
CV ← 0; Obj ← 0 19: for j ← 0 : N τ do 20: BO ← 1; TObj ← 0; X em ← X 0 21: for t loop ← 0 : h : T l do 22: Ξ ∼ N (0, Σ(X em )h) 23: X em ← X em + hf (X em , U i , Ξ, a x,0 ) 24: if (((X em , U i ) / ∈ Ω (X,U) ) & BO) then 25: CV ← CV + 1; BO ← 0 26: end if 27: TObj ← TObj + h(Γ 1 0 J acc l l + Γ 2 0 J roll l l ) 28
3. SIM function:

(a) The lines from 19 to 30 executes the MC simulation for the i th input -U i . The simulation consists of two for loops. The outer loop runs the MC simulation N τ

(N τ = γ
Ns ) times to maximize the eciency of CUDA cores to simulate several sce- nario instances. The inner loop is dedicated for problem (4.15), where the stochastic ODE is simulated by means of Euler-Maruyama integration and in due course of simulation, if the constraints are violated, the counter variable CV registers the violation and also, the objective is numerically approximated by means of Riemann sum. The inner loop objective is stored in a temporary variable TObj and the outer loop objective Obj is a N τ times accumulation of the inner loop objective. It is also important to note that the total number of MC simulations for each input is γ = N s × N τ and the total number of simulations for all the inputs is N = γ × n g .

(b) The line 31 returns the average objective and constraint violation certicate with respect to N τ simulations.

Scenario generation

As the SNMPC optimization problem in (4.15) is numerically solved by means of MC simulations, it is of paramount importance to sample enough number of scenarios to approximate the solution for the problem. Thanks to the theory of statistical learning and randomized algorithms which provides a systematic framework to derive the minimum number of scenarios required to achieve a probabilistic bound over the quality of the solution. The following content of this section is a summary of the core result presented in the seminal paper [START_REF] Vidyasagar | Randomized algorithms for robust controller synthesis using statistical learning theory[END_REF]]. To present the result more concretely, the context of the problem is laid down rmly with the following setup. Consider the task of approximating a stochastic function dened by

h(y) = E x∼Px [ψ(x, y)] (4.16)
where, x ∈ X is distributed w.r.t. the distribution P x over X and y ∈ Y. Let the empirical mean approximation of the function in (4.16) be dened with ĥ(y). Let the parameters for accuracy, probability level and condence level be dened with , β, δ ∈ [0, 1]. Given these parameters, the study of randomized algorithms is to derive a lower bound for the number of 4.5. Scenario-stochastic pNMPC scheme for control of semi-active suspension system 85 scenarios required to achieve the following requirement.

P[P[| ĥ(y) -h(y)| > ] ≤ β] ≥ 1 -δ, ∀y ∈ Y (4.17)
In simple words, the above (4.17) two layered probabilistic statement implies that the dierence between the empirical approximation and real function must be not greater than with a probability less than β and this must be certied with a probabilistic condence level of 1 -δ.

Given this prelude, the task of optimizing the function ĥ(y) by means of scenarios is dened by the following Theorem 4.5.1.

Theorem 4.5.1. ghoose the integers n nd m de(ned with

n ≥ ln( 2 δ ) ln( 1 1-β ) nd m ≥ 1 2 2 ln 4n δ (4.18)
qenerte iFiFdF smples y 1 , . . . , y n ∈ Y nd x 1 , . . . , x m ∈ X ording to P x F he(ne

ĥ(y i ) = 1 m m j=1
ψ(x j , y i ), i = 1 . . . n nd ĥ * = min 1≤i≤n ĥ(y i )

hen with on(dene 1 -δ it n e sid tht ĥ * is proly pproximte ner the minimum h(.) to ury nd level βF he result is universl nd pplile for ll fmily of funtions idysgr PHHIF

Using the results from Theorem 4.5.1, the empirical means in the SNMPC problem (4.15) both in the objective as well as the chance constraints can be numerically approximated.

An important point to note is that the chance constraint can be recasted with expectation formulation with

P[(X(.), U(.)) / ∈ Ω (X,U) ] = E[1{(X(.), U(.)) / ∈ Ω (X,U) }] (4.19)
where, 1{A} represents the indicator function over the set {A}. Thus, by leveraging the results propounded in Theorem 4.5.1, the number of scenarios γ for Algorithm. 2 can be derived. Setting β = δ = 0.05 yields the total number of input parameterization to n g = 64.

Utilizing the previous result and setting = 0.125, the number of scenarios (m = γ) for each parameterized input is γ ≈ 270. In total, the GPU simulates approximately 270 × 64 simulations over the prediction horizon at every sampling period. In this study is assumed to be xed, however it is really important to note that the parameter is implicitly a function of the computing resource such that the method is RT doable. As from Table 4.4, it clearly evident that with increase in the hardware conguration a considerable reduction in the mean and maximum CT is observed. However, the natural sampling period (T s ) of the platform is 5 ms and the best of embedded NVIDIA boards -Jetson Xavier hovers around 6 ms. Despite the boards doesn't compute the optimal input within T s , the scenario parameters , β, δ (mentioned in Section 4.5.4.2) can be tweaked to meet RT requirements, however the quality of the solution would be deteriorated. The histogram for the computational time for dierent NVIDIA embedded platforms is shown in Fig. 4.7.

4.5. Scenario-stochastic pNMPC scheme for control of semi-active suspension system 87 4.5.5.

Pareto optimality of objectives

As mentioned in Section 4.5.3, the two objectives (comfort and ride handling) are conicting in nature. Thus, the convex weights Γ 1 0 and Γ 2 0 ought to be tuned in order to strike a proper balance between the two objectives. To study the eects on variations of the convex weights, a Pareto optimality analysis was performed for several simulations. The simulation involved the vehicle moving at a constant velocity v x = 20m/s on a IS0-C road prole for a duration of 10s. In Fig. 4.8, the Pareto optimal front is plotted and this aids in detecting the right weights for the best results. The RMS values of the chassis acceleration for all the systems are listed in Table 4.5.

It is clearly evident that the RMS value of the proposed method is less than the nominal passive system and in par with the minimum passive system (Comfort design). For illustration purpose, Fig. 4.12 displays the chassis acceleration plot for ISO-C road prole. Also, at the same time, the proposed method minimizes the ride handling objective judiciously in comparison with other passive systems. For illustration purpose, Fig. .12: Chassis acceleration (z s ) plot for dierent controllers for ISO-C road prole of a massively multi-core processor for implementing the proposed parallelized version of the pNMPC controller. In this chapter, two methods has been proposed for control of vertical dynamics of a half-car vehicle through ER semi-active suspension system. The rst method proposes the parallelized pNMPC method where the road model is not accounted in the pNMPC problem formulation. In the second method, the stochastic road model in included into the pNMPC problem formulation and this method is termed the SS-pNMPC scheme.

The methods were both tested in simulation and as well as on NVIDIA embedded boards 90 Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active suspension system to validate and verify the feasibilty of the proposed approach. For future implementation of the parallelized pNMPC scheme under the purview of the EMPHYSIS project, the eFMU container could be potentially designed to be amenable on GPUs or any multi-core processors by the likes of GPUs. Thus, by this feature, several simulations of the eFMU models could be conducted in parallel and the parallelized pNMPC method could be seamlessly implemented for RT control.

Part III pNMPC -A code generation software tool for implementation of derivative free pNMPC scheme for embedded control systems 94 Chapter 5. pNMPC code generation tool 5.1 Introduction

Prelude

Over the last decade, there has been a tremendous amount of eort in the nonlinear MPC (NMPC) and optimal control community to build real-time (RT) operable software (S/W) for embedded control of engineering systems. There have been several contributions on code generation based NMPC/Optimal control toolboxes, however, by and large most of the existing toolboxes shares a common thread which is the need for derivatives of objective, constraints and dynamical model for solving the underlying optimization problem for the NMPC controller. In cases where the derivatives are not available and automatic dierentiation may not be feasible, the only recourse is to compute the derivatives using numerical methods such as nite dierence methods. However, numerically computing the derivatives is highly prone to error and as well as computationally taxing which could encumber RT feasibility. This brings the need for implementation of a complete derivative free NMPC package to address situations as mentioned before. As of today, the only derivative free NMPC toolbox available on market is the PDF-MPC package [START_REF] Alamir | tiliztion of nonliner systems using reedingEhorizon ontrol shemesX prmetrized pproh for fst systems[END_REF]] which uses MATLAB as the front end and MATLAB coder on top to deploy the code into embedded devices. The pNMPC code generation software tool presented in this chapter is highly inuenced from the aforementioned tool with improved support and features to cater the open source community and embedded control programmers. The pNMPC S/W is available in GitHub repository [START_REF] Rathai | pxwgX e gode qenertion ool por smplemenE ttion of pxwg gontroller por imedded gontrol ystems[END_REF]].

Motivation

Optimization solver plays a pivtol role in solving the underlying OCP for the MPC controller.

The entire edice of MPC controller reposes on the underpinning optimization problem. Typically, any state of the art optimization problem can be broadly classied into three types i) Zero order methods, ii) First order methods and iii) Second order methods. The order determines the type of information that the optimization solver can query to determine the solution, i.e. for zero order methods, only the function value can be obtained, for rst order methods, the function value and it's rst order derivatives (Jacobian) can be obtained and for second order methods, the function value, the rst (Jacobian) and second order derivatives (Hessian) can be obtained or approximated. In optimization parlance, this is demarcated by classifying zero-order methods as derivative free methods or black box optimization (BBO) problem and the rest as derivative based methods. The focus of this chapter circles around implementation of MPC controller using derivative free methods. Some of the compelling reasons for a need of a BBO based MPC are It is not uncommon in the real world industrial application, that the plant model and parameters are secured due to intellectual property (IP) rights and perhaps this sensitive information is never divulged even to the internal engineers working with the system.

Thus, the knowledge of the system is completely obscured from the control engineer and renders nearly impossible to design a MPC controller that predominantly relies upon the derivative information.

In cases where the model of system exists only as computer codes or in binary format (executable les), it becomes very cumbersome to obtain accurate derivatives by numerical methods. To exacerbate, if the code involves several break, if-else or goto statements, it is highly impractical to obtain the derivatives using automatic dierentiation methods. Case in point, this is common in the functional mock-up interface (FMI) standard [START_REF] Blockwitz | Functional mockup interface 2.0: The standard for tool independent exchange of simulation models[END_REF], which is used for model exchange and co-simulation purposes.

Also, this is the main framework of the EMPHYSIS project [EMPHYSIS 2017].

Even when the functional form is available to the control engineer, at times, computation of derivatives can be computationally too expensive or noisy and this could preclude from practical implementation of MPC controller for fast sampled systems. Also, in cases where the function is discontinuous and not dierentiable, the derivative based methods can lead to undened behavior.

In today's world, with availability of deluge of data and increased computation power, data-driven modeling has challenged the pre-existing notions of rst principles modeling and perhaps it wouldn't be too much of a stretch to consider the former superseding the latter in a few decades now. Machine learning models such as neural networks, Gaussian processes etc. provide excellent empirical approximations of the underlying dynamical systems and typically one can right away utilize these models for control design. It is also important to note that many of these models are not dierantiable by denition.

As per the cited reasons, it is certainly not an unreasonable requirement but, also, a matter of paramount importance to address the control problem by designing a derivative free MPC scheme. In this chapter, a BBO based parameterized NMPC (pNMPC) S/W tool is proposed to circumvent the aforementioned issues. The term parameterized refers to the parameterization of the control input [START_REF] Alamir | tiliztion of nonliner systems using reedingEhorizon ontrol shemesX prmetrized pproh for fst systems[END_REF]]. The parameterization serves two purposes a) reduces the computational burden for the optimization solver and b) design of a parsimonious control input prole. It is also important to note that the potential of the pNMPC method as well as the toolbox relies upon the astuteness of the control engineer to model the input prole with ecient, eective and economical parameterization.

Related works

There has been several contributions on development of code generation based MPC toolbox for embedded systems. The rst automatic NMPC code generation toolbox traces back to

AutoGen [START_REF] Ohtsuka | Automatic code generation system for nonlinear receding horizon control[END_REF] The µAO-MPC [START_REF] Zometa | µAO-MPC: a free code generation tool for embedded real-time linear model predictive control[END_REF] toolbox provides code generation feature for embedded RT linear MPC. The underlying optimization solver for the tool involves a rst order method based on Nesterov's gradient method coupled with augmented Lagrangian method. The toolbox provides front end interface to Python and the generated C code is implemented on embedded systems.

The CVXGEN [START_REF] Mattingley | CVXGEN: A code generator for embedded convex optimization[END_REF]] toolbox provides code generation feature for an embedded convex optimization problem. The front end involves a quadratic programming (QP) modeling framework, where the user enters the optimization variables, dimensions, parameters, objective, inequality/equality constraints etc. and the output of the program is an optimized C-code of the optimization problem. The toolbox exploits the structure of the QP problem such as sparsity, optimal KarushKuhnTucker (KKT) matrix factorization etc. to speed up the computation. Despite the toolbox is agnostic to any specic application, the method has been mostly used for MPC controller design. The only limitation is that the method applies only for a convex QP problem.

FiOrdOs [START_REF] Ullmann | FiOrdOs: A Matlab toolbox for C-code generation for rst order methods[END_REF]] is an automated C-code generation MATLAB toolbox for rst order methods for parameteric convex programs. The only limitation of the toolbox is that it can only handle only convex QPs with simple convex sets on which projections can be evaluated at low cost.

The ACADO toolbox [START_REF] Houska | An auto-generated real-time iteration algorithm for nonlinear MPC in the microsecond range[END_REF] provides an automatic C-code generation feature for implementation of NMPC controller based on real-time iteration (RTI) using Gauss-Newton method. The toolbox is implemented in C++ and also provides interface with MATLAB. The front end involves a modeling framework to model the states, inputs, parameters, dierential equations, equality/inequality constraints, integrator parameters, solver parameters etc. and the back end involves the optimization solver, which is outsourced to other third party S/W such as qpOASES [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], QPdunes [START_REF] Kouzoupis | Block condensing for fast nonlinear MPC with the dual Newton strategy[END_REF], FORCES [START_REF] Zanelli | FORCES NLP: an ecient implementation of interior-point methods for multistage nonlinear nonconvex programs[END_REF], IPOPT [Wächter and Biegler 2006b]. The nal output is the respective C codes of the model, constraints, integrator sensitivities etc.

The VIATOC toolbox [START_REF] Lozoya-Santos | A toolkit for nonlinear model predictive control using gradient projection and code generation[END_REF] which shares similar syntax as ACADO toolbox, however the NMPC problem is specically tailored to an optimization solver based on projected gradient method. The toolbox also provides automatic code generation feature similar to ACADO. The toolbox was programmed in C++ environment.

The Multi-Parameteric toolbox (MPT) [START_REF] Kvasnica | Automatic code generation for real-time implementation of model predictive control[END_REF] provides code generation feature for implementation of explicit MPC. The method uses QP parameteric programming method and precomputes the optimal control input over the partitioned polyhedral regions of the constraint space. The states are considered as the parameters and the optimal feedback law given current state of the system is nally computed by performing a table traversal algorithm using binary tree search method and the optimal control input is recovered. The method is suitable for linear MPC problems with low state, input dimensions and horizon length. The toolbox is programmed in MATLAB and uses MATLAB coder to generate C-code for embedded systems.

The SPLIT toolbox [START_REF] Shukla | Software and hardware code generation for predictive control using splitting methods[END_REF]], a C code generation tool for linear MPC problems uses operator splitting methods and also, the toolbox is capable of generating both software as well as hardware specic code such as FPGA boards. The toolbox provides a MATLAB front-end interface.

The Optimization Engine toolbox (OpEn) [START_REF] Sopasakis | OpEn: Code Generation for Embedded Nonconvex Optimization[END_REF]] is a code generation tool for RT embedded non-convex optimization problems. The OpEn tool combines proximal averaged Newton-type method with penalty and augmented Lagrangian methods.

The method was implemented in Rust programming language and provides front end interfaces with MATLAB, Python, C/C++ or via a TCP socket.

The ParNMPC [START_REF] Deng | A parallel code generation toolkit for nonlinear model predictive control[END_REF]] is a parallel code generation toolbox to generate C/C++ code for NMPC controller. ParNMPC utilizes the OpenMP programming framework and implements Newton-type methods across the multi-core processors.

Chapter contributions

The main contribution of this chapter is to present a derivative free pNMPC toolbox for embedded control of engineering systems. The salient features of the pNMPC S/W are

The pNMPC S/W was completely programmed in C++ environment. The S/W by it's inbuilt design provides an OCP modeling framework, where the user can provide the objectives, constraints, dierential equations, Black box modules, solver parameters, OCP parameters in a hassle-free way.

The S/W is completely derivative free and as well as free of any online matrix operations such as matrix-matrix addition, multiplication, inverse computation, linear algebra solve etc.

A one-stop control solution is provided to the end user for generation of portable, efcient, optimized and embeddable C code of the pNMPC controller. Extensions to MATLAB/Simulink environment is also provided. The inbuilt code generation module is scrupulously engineered, thus leading to more optimized code in terms of memory footprint and computation time.

The toolbox is independent of any external libraries such as Basic Linear Algebra Subprograms (BLAS), Linear Algebra Package (LAPACK) etc. and consumes less memory footprint.

The S/W also includes Graphic Processing Unit (GPU) based parallel implementation of the optimization solver as well as Compute Unied Device Architecture (CUDA) code generation module suited for problems with large optimization variables.

Additional to the inbuilt modeling framework, the S/W also provides the exibility to include and call external functions (black box functions) which can be in the form of either source code or static/dynamic libraries.

The S/W is open sourced under GNU-LGPL v3 license.

The proposed pNMPC S/W tool was tested in Hardware in the Loop (HiL) simulation for a quarter car vertical dynamics model on dSPACE MicroAutoBoX-II (MABXII, See Section 5.7).

5.2 pNMPC theoretical background

pNMPC problem formulation

The goal of the pNMPC toolbox is to solve the following Optimal Control Problem (OCP) problem at every sampling period. min x(.), p(.)

T 0 l(x, u(x, p, κ), κ) dt + ψ(x(T ), κ(T )) subject to ẋ = f (x, u(x, p, κ), κ, t), ∀t ∈ [0, T ] u(x, p, κ) ∈ U, x ∈ X , p ∈ P , ∀t ∈ [0, T ] x(0) = x 0 , x(T ) ∈ X T (5.1)
where, x ∈ R nx , u ∈ R nu , p ∈ R np and κ ∈ R nκ represents the state vector, input vector, input parameterization vector and external parameters (i.e. model parameters, set point for tracking or measured disturbances) vector respectively. The input map u : R nx ×R np ×R nκ → R nu maps the states, input parameterization vector and external parameters to actual input for the system. The sets X , U, P and X T denote the state constraint, input constraint, input parameterization constraint and the terminal state constraint respectively. The OCP is subjected to a set of dierential equations denoted with f : R nx × R np × R nκ × R + → R nx and the Lagrangian cost (stage cost) and Mayer's cost (terminal cost) are denoted with l :

R nx × R np × R nκ → R and ψ : R nx × R np × R nκ → R respectively.
It is important to note that the cost terms can be economic and no stipulations are enforced such as non-negativity or convexity. It is also important to note that handling the impact of the objectives or constraints design or choise of sampling period (∆t) on the closed loop behavior is left to the user.

The pNMPC S/W transcribes the above OCP formulation (5.1) into a generic constrained optimization problem. The OCP transcription is implemented by discretizing the problem in time with a nite time step of ∆t. Thus, the states of the system are eliminated from the problem by simulating the dierential equation over time (direct single shooting method, see

Chapter 2), the integral objective is numerically approximated (Riemann sum) and the constraints are quantied at every discretized time step. Let J and h i , ∀i ∈ {1 . . . n c } denote the transcribed objective and inequality constraints (n c inequality constraints). The transcribed constrained optimization problem is described by the following form where p ∈ R np is the optimization vector or decision variables which is also the input parameterization vector. Let the solution of the above problem (5.2) be denoted with p * . Utilizing the optimized input parameterization vector, the optimal input u(x * (τ ), p * (τ ), κ(τ )) is injected into the system over the time period τ ∈ [0, ∆t]. Henceforth, by marching forward in time and with receipt of new state vector, this process is repeated in receding horizon manner.

Visualization of control parameterization

To understand the representative and expressive power of the input parameterization technique, consider a sinusoidal input parameterization over a prediction horizon T = 2s as described below.

u(p(t), t) = p 1 (t)sin(2πp 2 (t)t + p 3 (t)) (5.3) where the time dependent input parameterization vector is p(t) p(t) ∈ P. The physical interpretation of the input parmeterization vector corresponds to the amplitude, frequency and phase of the sinusoid input respectively. In order to discretize the continuous time input parameterization vector p(t), consider n p 1 , n p 2 , n p 3 nite control points over the prediction horizon, then each of parameterized input at a specic time instant can be described with the pairs {t i , p 1,i }

i=np 1 -1 i=0 , {t i , p 2,i } i=np 2 -1 i=0 , {t i , p 3,i } i=np 3 -1 i=0
. For a simplistic case, let the control points be connected with a linear interpolation curve and let n p 1 = 5, n p 2 = 6, n p 3 = 8, then the respective input parameterization proles and the sinusoid control input is illustrated in Fig. 5.1 and Fig. 5.2.

By virtue of these parameterization points, a exible control input prole can be designed and in the context of pNMPC toolbox, the parameterization control points are the optimization variables which ought to be ascertained from the solver. In this sinusoid example, in total there are n p 1 + n p 2 + n p 3 = 19 control points, which in this example is certainly an over parameterization of the control input. However, in practice, with a careful and deliberate placements of control points one can obviate the perils of over-parameterization of control input.

Derivative free optimization module

In this section, the underlying BBO module for the proposed pNMPC controller is discussed in detail. The content and core results summarized in this section is based on the Sequential Quadratic Programming based Black Box Optimization (SQP-BBO) method proposed in [START_REF] Alamir | tiliztion of nonliner systems using reedingEhorizon ontrol shemesX prmetrized pproh for fst systems[END_REF]], [START_REF] Alamir | tiliztion of nonliner systems using reedingEhorizon ontrol shemesX prmetrized pproh for fst systems[END_REF]], [START_REF] Alamir | tiliztion of nonliner systems using reedingEhorizon ontrol shemesX prmetrized pproh for fst systems[END_REF]]. The method falls under the category of interpolation based trust region methods for derivative free optimization. The method is based on the technique of sequentially approximating the cost and constraint functions by means 5.3. Derivative free optimization module 101 of quadratic functions and at successive iterations the optimal solution is computed based on multiple trust region switch conditions. In the Subsection 5.3.2, uni-variate case of the optimization problem is discussed in detail and in Subsection 5.3.3, the method is extended for the multi-variate case.

Constraint reformulation (Scalarization)

Consider an optimization problem as dened below min p∈R np J(p) s.t. g(p) ≤ 0 (5.4) where p ∈ R np is the vector of optimization variables and J and g represents the scalar cost and scalar constraint of the optimization problem respectively. In cases where there exists several inequality constraints, then all the constraints are scalarized by either of the two following forms.

Consider there exists n c constraints acting upon the optimization problem, i.e. where the optimization variable p ∈ R belongs to a bounded interval [p min , p max ] with p max ≥ p min . In order to dene a local quadratic approximation of a function f (f is a generic representation of a function which can be either J or g) over an interval I, consider a variable α > 0 with respect to a point p such that the interval I is dened with The above points are equivalently represented with {f -, f 0 , f + } respectively. The locally approximated quadratic function q f (p) can be expressed using a parabolic parameteric form with {f -, f 0 , f + }, the coecients can be computed by solving the following linear algebra problem.

I := [p -α, p + α] ∩ [p min , p max ]
q f (p) = a f p -p c β 2 + b f p -p c β + c f
  1 -1 1 0 0 1 1 1 1     a f b f c f   =   f - f 0 f +   (5.13)
The solution for the above problem (5.13) yields the results

a f = 1 2 [f -+ f + ] -f 0 b f = 1 2 [f + -f -] c f = f 0 (5.14)
In order to compute the local minimizer for the approximated quadratic function ( 5.12), it is of high importance to consider the parameter a f to determine the existence of minimizer, i.e. when a f = 0, nite value for the minimizer exists and when a f = 0, the solution can be either at -∞ (when sign(b f ) > 0) or +∞ (when sign(b f ) < 0). Thus, the solution p

(f ) s is expressed with p (f ) s =          p c - βb f 2a f , if a f = 0 +∞, if a f = 0 and b f ≤ 0 -∞, if a f = 0 and b f > 0 (5.15)
In order to conne the solution (5.15) obtained from the local quadratic approximation within the interval I dened in (5.10), the solution is projected onto the interval I and the projected solution is dened with

p (f, * ) s := min p max , max p min , p (f ) s (5.16)
The optimal function value of the locally approximated quadratic function at the projected solution (5.16) is dened with

q * f = a f p (f, * ) s -p c β 2 + b f p (f, * ) s -p c β + c f (5.17)
Using the computed q * f value from (5.17) the extreme values of parabola q f (.) over the interval of interest I are obtained. These extreme values are dened with 104 Chapter 5. pNMPC code generation tool

q min f := min f -, f + , q * f q max f := max f -, f + , q * f (5.18)
The values of p at these extreme values are denoted with p min

f :=        p min , if q min f = f - p max , if q min f = f + p (f, * ) s , if q min f = q * f p max f :=        p min , if q max f = f - p max , if q max f = f + p (f, * ) s , if q max f = q * f (5.19)
When the generic function (f ) is the inequality constraint function i.e. f = g, special cases arises for the local quadratic function (q g (p)) approximation over the interval I. The possible cases are 1. q g (p) is non-negative for all p ∈ I and this applies when q min g > 0. In this case, g is non-negative in the interval I This means that the inequality constraint is violated over the entire interval I.

2. q g (p) is negative for all p ∈ I and this applies when q max g ≤ 0. This means that the inequality constraint is admissible over the entire interval I.

3. q g (p) is negative on a strict subset of the interval I and this occurs when q min g ×q max g < 0.

In the last case, there could be either one or two values of p that belong to the interval I and this occurs only when q g (p) = 0 and these value can be obtained by analyzing the discriminant ∆ := b 2 g -4a g c g and when ∆ ≥ 0 (non-negative) there is at least one real solution. Thus, the possible candidate solutions are p (0,+)

g := p c + βmax -b g ± √ ∆ 2a g p (0,-) g := p c + βmin -b g ± √ ∆ 2a g (5.20)
Let Z - g ⊂ I denote the subset of values of p that belong to I where q g (.) is negative. The set Z g is computed with 5.3. Derivative free optimization module 105

Z - g :=                      [p min , p (0,-) g ], if p (0,-) g < p min & g -≤ 0 [p (0,+) g , p max ], if p (0,-) g < p min & g + ≤ 0 [p min , p (0,-) g ], if p (0,+) g > p max & g -≤ 0 [p (0,-) g , p max ], if p (0,+) g > p max & g + ≤ 0 [p (0,-) g , p (0,+) g ], if [p (0,-) g , p (0,+) g ] ⊂ I & a g > 0 A ∪ B, if [p (0,-) g , p (0,+) g ] ⊂ I & a g < 0 (5.21)
where, A = [p min , p (0,-) g

] and B = [p

(0,+) g , p max ].
It is important to note that the set Z g is either an interval or a union of two intervals and this dierence is demarcated with n z g ∈ {1, 2} which corresponds to an interval I

(1) g or union of intervals I

(2) g . All the terminologies used are summarized and tabulated in Table 5.1.

p c Center I = [p -α, p + α] ∩ [p min , p max ] β Semi-length of I q f (.) Local quadratic approximation of f over I a f , b f , c f Coecients of parabola q f (.), f ∈ {J, g} p (f ) s Position of singular point q f (.) p (f, * ) s Projection of p (f ) s over I q * f
The value of parabola q f (.) at p (f, * ) s q min f Minimum value of q f (.) on I q max f Maximum value of q f (.) on I p min f Location of minimum value of q f (.) on I p max f Location of maximum value of q f (.) on I p

(0,+) g , p (0,-) g Solution of q g (p) = 0 Z - g Subset of I where q g (.) ≤ 0 n z g Number of intervals in Z - g I (1) g , I (2) g Interval dening Z - g
The BBO algorithm is completely premised upon the previously dened three cases. The algorithm for the next iteration p (i+1) and trust region update size α (i+1) are given by the following steps 1. When q max g ≤ 0, then the whole interval I is the search space and the function J is minimized over the whole interval I. The candidate value of the update p (i+1) is given by p cand ← p min J (5.22) The minimizer in the above equation (5.22) is obtained from equation (5.19) where the function f is replaced with J. This computation is based on the assumption that the 106 Chapter 5. pNMPC code generation tool quadratic approximation is appropriate. The logical condition to verify this assumption and also, to update the trust region size is given by C ← C 1 C 2 where,

C 1 ← J(p cand ) < J(p (i) ) g(p cand ) ≤ 0 C 2 ← J(p cand ) ≤ J(p (i) ) g(p cand ) < 0 (5.23)
2. When q max g > 0, which means that the constraints are strictly non-negative which tantamount to constraint violation and the priority ought to be given to minimization of the inequality constraint g. In this case, the candidate value of the update p (i+1) is given by p cand ← p min g (5.24)

and the trust region update condition is given with C ← g(p cand ) < g(p (i) )

(5.25)

3. When q max g ≥ 0 and q min g ≤ 0, which means that there exists a subset in I where there exists a solution. In this case, the integer n z g and the corresponding intervals I

(1) g and I

g are computed. Utilizing these intervals, the potential candidate for p (i+1) update is computed by

p cand ← arg min p∈{p (min,l) J } n z g l=1
J(p) (5.26) where p (min,1) J and p (min,2) J are optimal solutions that minimize J over the intervals I

(1) g and I

(2) g respectively. The trust region update condition is given by

C ←    g(p cand ) < g(p (i) ) , if g(p (i) ) > 0 J(p cand ) < J(p (i) )
g(p cand ) ≤ 0, else (5.27) The update of next iterate p (i+1) and α (i+1) is implemented according to following rules.

If C is true, then p (i+1) is assigned to the computed candidate value p cand . The trust region parameter α (i) is increased according to

α (i+1) ← β + • α (i) ; β + > 1 (5.28)
Otherwise, the current value p (i+1) ← p (i) is used and the trust region size is decreased according to

α (i+1) ← max α min , β -• α (i) ; 0 < β -< 1 (5.29)
Let N iter represent the number of iterations the above algorithm is repeated, then the total number of function evaluations of the objective and constraint functions (J, g) is N eval = 4N iter + 1 The multi-variate case is an extension to the uni-variate case, where an uni-variate optimization problem is solved over each component of the decision variables while the rest of the values are maintained constant. Consider the decision variables to be p ∈ R np and let the list of decision variables be indexed with l ∈ {1, 2, . . . , n p }. Let η ∈ R denote the scalar variable over which the uni-variate optimization is performed. In notational form, this is expressed with p (η,l) and an element in R np is dened as p (η,l)

j := p j if j = l η if j = l (5.31)
The formulation (5.31) is extended for ∀j ∈ {1, 2, . . . , n p }. The total number of loop count to visit all the n p components for N iter iterations and N eval is given by

N loop = N eval -1 4n p × N iter (5.32)
It is important to note that a feasible choice of pair (N eval , N iter ) and must satisfy the inequality N eval ≥ 4n p N iter + 1 OCP design and specication: The OCP's design and specication are programmed by the user using C++ as a front-end modeling language.

Lexical analyzer: Lexical analyzer or tokenizer takes the user's OCP design and specication and breaks down the entries into separate characters or special tokens such as the states, inputs, parameters, math operations etc.

Syntax analyzer: Syntax analyzer or parser scans through these tokens and contrives a relational tree or parse tree. Consider an example y = x 1 x 2 + sin(x 1 x 2 ), then the computed parse tree for this relation is illustrated in Fig. 5.5.

x 1 x 2 x 1 x 2 * * Sin + Code optimization: The code optimization stage is a crucial stage of the code generation module. Typically, the code optimization module obliterates redundant relations, self negation operations, identity and inverse relations for addition, subtraction, multiplication and division operations such as (to list a few) 0 * x, x + 0 etc. In the above example, it is clearly evident that the term x1 * x2 is computed twice. After code optimization, this redundancy is removed and the stack is updated.

Target C code generation: In this stage, the embedded C les are exported by using the le stream operations. Under circumstances of securing the generated source code, either a static or dynamic library can be created. However, this has to be done manually by the user.

Simulink S-function: This stage is optional, however in today's world, embedded control has virtually become MATLAB/Simulink's efdom and it plays a dominant role in design, development and deployment of production code to embedded systems both in industry and academia. The pNMPC S/W provides the exibility to provide Simulink compliant C codes and this can be included into Simulink by availing MEX wrappers in MATLAB or C-MEX S-functions features from Simulink.

Embedded system: Finally, the generated code is deployed into the embedded system either through Simulink or manually by the user.

Application of pNMPC toolbox

The pNMPC S/W was tested for several examples and the simulation results looks promising and viable for RT implementation. In this Chapter, in the interest of space, the results and simulation of two examples are described in detail which are 1. Cart-pole swing up problem.

PVTOL stabilization problem.

The outputs were compared against ACADO toolkit [START_REF] Houska | An auto-generated real-time iteration algorithm for nonlinear MPC in the microsecond range[END_REF]

Cart-pole swing up problem

The task of the cart-pole swing up problem [START_REF] Mills | Nonlinear model predictive control of an inverted pendulum[END_REF] is to stabilize a pole in an upright direction (typically starting from downward position) which is attached to a movable cart by means of a revolute joint. The control input to this system is the horizontal force applied to the cart and system is bounded by physical constraints which are the length of cart travel and the input force. The nonlinear state space equations of the system are given below

ẋ = v v = -m 2 lsin(θ) θ2 + u + m 2 gcos(θ)sin(θ) m 1 + m 2 (1 -cos 2 (θ)) θ = ω ω = -m 2 lcos(θ)sin(θ) θ2 + ucos(θ) + (m 1 + m 2 )gsin(θ) l(m 1 + m 2 (1 -cos 2 (θ))) (5.34)
where m 1 , m 2 , l, g represents the mass of the cart, mass of the pole, length of the pole and acceleration due to gravity respectively. The values for the parameters are listed in Appendix A. The state vector of the system are x = [x, v, θ, ω], which are the cart position, cart velocity, pole angle and pole angular velocity respectively and input to the system is u, which represents the force acting on the cart. The constraints acting on the system are 5.35) where x max and u max represents the maximum bounds of the cart position and cart force respectively. The OCP for the system is given as min x(.), p(.) 5.36) where t f , Q f , Q, R and T s represents the look ahead period, quadratic terminal state cost, quadratic stage state cost, quadratic input cost and sampling period respectively. The input parameterization is an ane state feedback policy where the parameterization vector is p = [p 1 , p 2 , p 3 , p 4 , p 5 ] which is modeled a constant over the horizon. Once the OCP is solved for p * , the input u(p * , x(0)) is injected into the system over the period T s and this process is repeated in a receding horizon fashion. To compare against the ACADO controller, an un-parameterized version of the NMPC problem (5.36) was implemented with the following settings -Integrator -4th order Implicit Runge Kutta integrator, QP solver -qpOASES, Hessian approximation -Gauss-Newton, Discretization -Multiple shooting, Discretization intervals -30 and for the rest, default parameters were utilized. The study was conducted in two parts with two initial conditions which are x(0) = [0, 0, π 2 , 0] and x(0) = [0, 0, π, 0]. pNMPC and ACADO controller for the second case respectively. From the plots, it is evident that the ACADO controller crashes, however, despite the numerical ill-conditioning of model, the pNMPC controller fares well with state feedback parameterization and at the same time, the system is stabilized. Fig. 5.9 illustrates the parameterization values for the second case.

-x max ≤ x ≤ x max -u max ≤ u ≤ u max (
x(t f ) T Q f x(t f ) + t f 0 (x T Qx + u T Ru) dt subject to (5.34), (5.35), x(0) = {[0, 0, π, 0], [0, 0, π 2 , 0]} u(p, x) = p 1 x + p 2 v + p 3 θ + p 4 ω + p 5 (

PVTOL stabilization problem with Black-box models

The task of the PVTOL (planar vertical takeo and landing aircraft) [START_REF] Martin | A dierent look at output tracking: control of a VTOL aircraft[END_REF] stabilization problem is to regulate the states of the system to the origin given where, the state vector is x = [y, v y , z, v z , θ, θ] which represents the vertical position, vertical velocity, horizontal position, horizontal velocity, roll angle and roll rate respectively, the input vector is u = [u 1 , u 2 ] which represents the lift acceleration and angular acceleration respectively and the parameter vector is κ = [σ] which represents the coupling between roll and lift eects. The variables ψ BB y (x, u, κ) = -u 1 sin(θ) + σu 2 cos(θ) and ψ BB z (x, u, κ) = u 1 cos(θ)+σu 2 sin(θ)-1 represents the dynamics of the system which are deliberately modeled as a black-box model and invoked using the function calls "ModelY" and "ModelZ" respectively.

The input arguments for these functions are the state vector, input vector and parameter vector and the output is the respective dynamics of system. The source codes for these functions were compiled to a dynamic library le and linked during the compilation process of the pNMPC controller. This example serves as an use case of the pNMPC S/W, where external black box models can be linked with the S/W's inbuilt symbolic variables. , θ max represents the maximum bounds over the inputs u 1 and u 2 and the roll angle state θ. The OCP for the stabilization problem is dened as min x(.), p(.) (5.37), (5.38), x(0) = [1, -1, 2, 1, -14π 10 , -0.1] u(p) = p (5.39) where t f , Q, R and T s represents the look ahead period, quadratic stage state cost, input cost and sampling period respectively. The input parameterization is p = [p 1 , p 2 ] where each parameter has two control points placed equidistantly over the prediction horizon and follows a linear prole. Once the OCP is solved for p * , the input u(p * ) is injected into the system over the period T s and this process is repeated in a receding horizon fashion. The following setting was chosen for the ACADO controller, Integrator -4th order Runge Kutta integrator, QP The constraints imposed on the system were the min/max damper force constraints, the min/max stroke deection constraints and the min/max PWM-DC (Pulse Width Modulation -Duty Cycle) signal which operates the damper characteristics.

x(t f ) T Q f x(t f ) + t f 0 (x T Qx + u T Ru) dt subject to
The objective was considered to maximize the comfort, which in turn is reected in minimizing the chassis acceleration of the vehicle. The method deserves some explanation to elucidate it's working principle. Consider an example as illustrated in Fig. 5.16. The optimization vector is four and black dots represents the unoptimized values of the variables in its respective places at the current iteration. In the rst iteration, the optimization vector is disseminated to four CUDA processors and an univariate SQP-BBO routine is executed in parallel for each component of the optimization vector with rest of the components held constant with previous values. The cross mark represents the component of interest for optimization and as well as the optimized component. Once the optimization is carried out in the CUDA core, the optimized components are collected and the optimization vector is updated and this is passed on for the second iteration. This procedure is repeated till the end of iteration count. where, the state vector is x = [r, r d , θ, ω] which represents the crane position, velocity, crane angle of swing and angular rate of swing respectively. The input is u which represents the applied input force. M , m, l, f r and g represents the mass of the moving body, mass of the suspended mass, length of the crane, drag resistance and acceleration due to gravity respectively. The constraints acting on the system are -u max ≤ u ≤ u max -θ max ≤ θ ≤ θ max -ω max ≤ ω ≤ ω max folds lower than the pNMPC-CPU version and this dierence would be more pronounced when large parameterization vectors (greater than 128) are used, as the throughput of CUDA GPU 5.9. Future works and conclusions 127 devices are much larger. It is also important to note that the sampling period T s considered for this example is 0.1s and the mean computation time of pNMPC-CPU controller hovers around 0.2232s and renders impossible for RT implementation. However, the mean computation time of pNMPC-GPU controller hovers around 0.036s, which is feasible for RT implementation. Smart parameterization -As mentioned in the Section 5.1.2 and to re-emphasize the fact again that the potential of the method solely depends upon the parameterization technique adopted by the control engineer. However, in cases when this becomes intricate in nature, the onus of determining an optimal parameterization is substantially increased.

In such a situation, it would be remiss of not utilizing tools from the machine learning community. Methods developed in reinforcement learning (RL) [START_REF] Sutton | einforement lerningX en introdution[END_REF] zones in parallel to the proposed approach and by availing tools and methods from RL community a smart parameterization technique for pNMPC controller can be developed.

Implicit solvers -As of now, the proposed S/W provides only explicit ODE solvers for the pNMPC controller. However, in many applications where the system is intrinsically sti in nature, one is obliged to use implicit solvers for numerical stability. The future version of the S/W would encompass support for several implicit solvers.

Equality constraints -In the future, the S/W would include support for equality constraints, which is important for control of periodic systems or to enforce time point constraints on the system. pNMPC =>s e t I n i t i a l T i m e ( 0 ) ; pNMPC =>s e t F i n a l T i m e ( 2 ) ; pNMPC =>s e t S t e p S i z e ( 0 . 1 ) ; pNMPC =>g e t S o l v e r ()=> s e t N i t e r ( 4

N

  s vs J norm CL for dierent values of computational scale factor γ . . . . . . . . . 3.8 Control scheme for the proposed analysis . . . . . . . . . . . . . . . . . . . . . . 3.9 dSPACE MicroAutoBox II -1401/1511 . . . . . . . . . . . . . . . . . . . . . . . 3.10 Computational eciency between the pNMPC and linearization based MPC controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Chirp Road prole with amplitude of 1 mm and frequency sweep from 5 to 25 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 0.2 Introduction and problem statement of the thesis . . . . . . . . . . . . 1 0.3 Structure of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 0.1 Thesis framework This thesis is a compilation of the results obtained over three years of doctoral work (beginning from September 2017 to November 2020) on Embedded modeling and control of vehicle dynamics: application to a small car pilot plant with ER dampers conducted in the team SAFE (Safe, Controlled and Monitored Systems), formerly known as SLR (Systèmes Linéaires et Robustesse), Automatic department, GIPSA-lab, Grenoble. The thesis was conducted under the guidance of Prof. Mr. Olivier SENAME (Director of the thesis, Professor Grenoble INP) and Prof. Mr. Mazen ALAMIR (Co-director of the thesis, Research Director CNRS). This work was supported and funded by the ITEA3 European project, 15016 EMPHYSIS (EMbedded systems with PHYSical models In the production code Software) from 2017 to 2020 [EMPHYSIS 2017].

  Figure 1.1: CAD design of double wishbone vehicle suspension system

Figure 1 . 2 :

 12 Figure 1.2: Force-deection velocity characteristics (SER) for passive, semi-active and active suspension systems

Figure 1 . 3 :

 13 Figure 1.3: Quarter car vehicle model

  INOVE test platform is a GIPSA-lab experimental setup built in collaboration with SOBEN company and the platform is a 1:5 scaled 4 poster testbed which is exclusively dedicated to study and assess the performance of vertical dynamics of vehicle under several road prole scenarios. The test platform consists of three major components which are a) Host PC, b) Target system and c) Physical plant (illustrated in Fig. 1.5). In systems point of view, the test platform is used to verify and validate controllers, observers and also, perform fault diagnosis and condition monitoring [Vivas-Lopez et al. 2014a] of the components. The INOVE platform is shown in Fig. 1.6. 1.5.1 Hardware description Host PC -The host PC is where the controller/observer/fault diagnosis/condition monitoring modules are designed, developed, deployed and tested to verify and validate the 1.5. INOVE test platform (GIPSA lab) module's performance on the target system. The front-end S/W is MATLAB/Simulink and by virtue of Simulink Real-Time Workshop (RTW), the Simulink model is converted to embedded C code and deployed on the target. Target system -There are two target systems for real-time (RT) implementation which are xPC target and dSPACE MABXII. The targets run Simulink RT operating system. Data acquisition system (DAQ) -The I/O cards include 2x National Instruments PCI-6259 DAQ boards, 1x Quatex QSC 100 serial board. The sampling frequency of the DAQ is 200 Hz, i.e. a sampling period of 5 ms.

Figure 1 . 5 :

 15 Figure 1.5: Schematic of INOVE test platform [Vivas-Lopez et al. 2014a]

Figure 1

 1 Figure 1.6: INOVE test platform GIPSA-lab

Figure 1 . 7 :

 17 Figure1.7: Schematic of data ow between the simulation tool and the FMU[START_REF] Blockwitz | Functional mockup interface 2.0: The standard for tool independent exchange of simulation models[END_REF]

  Figure 1.8: Simplied schematic of eFMI workow [EMPHYSIS 2017]

Figure 1 . 9 :

 19 Figure 1.9: Detailed schematic of eFMI workow [EMPHYSIS 2017]

Figure 1 .

 1 Figure 1.10: Detailed schematic of eFMI workow [EMPHYSIS 2017]

[

  [START_REF] Nocedal | xumeril optimiztion[END_REF]. The most commonly used second order derivative based methods for solving NLP problems are a) Sequential Quadratic Programming (SQP) methods and b) Interior-Point (IPM) methods. The SQP method can be further classied based on its working principle which are a) SQP Active Set (AS) methods and b) SQP Interior Point Quadratic

  represents the dynamics of the system (Ordinary Dierential Equations (ODEs)/Dierential Algebraic Equations (DAEs)), inequality and equality constraints respectively. x 0 and r : R n → R r represents the initial and nal state conditions (terminal set) on the state variable. The inequality constraints

  2.2) for solving the non-linear OCP problem (2.5), which are Hamilton-Jacobi-Bellman (HJB) equation or Dynamic programming (DP) Direct methods or Direct transcription methods Pontryagin minimum/maximum principle (PMP)

Figure 2 . 2 :

 22 Figure 2.2: Classication of optimal control problems

Figure 2 . 3 :Fig. 2 . 3 .

 2323 Figure 2.3: Direct single shooting illustration

Figure 2 . 4 :

 24 Figure 2.4: Direct multiple shooting illustration

  benets of utilizing direct multiple shooting methods include a) Better simulator stability with unstable system, b) Parallelizability of ODE/DAE simulations, c) Structural properties of the Hessian matrices aid the optimization routine. The Jacobians and Hessians for the ODE/DAE simulator are obtained via either forward or adjoint sensitivity analysis methods

Figure 3

 3 Figure 3.1: Measured damper force (u) for PRBS based road prole

Figure 3 . 2 :

 32 Figure 3.2: PWM-DC signal (φ)

Fig. 3 .Figure 3 . 3 :

 333 Figure 3.3: Wavelet analysis of ER semi-active damper force (u) signal

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: F ER vs z def plot for dierent PWM-DC signals

Figure 3 . 6 :

 36 Figure 3.6: Predicted damper force and measured damper force u

  For the pNMPC design considered, six constraints are included in the problem formulation which are Semi-active ER damper input constraints: PWM input constraint: φ(t) ∈ U. Max/Min damper force constraint: This forms a non-linear mixed state-input constraint such that u(t) ∈ [u, u], where u and u are the minimum and maximum saturation forces for the semi-active suspension system. State limitations constraints: Max/Min deection between the chassis and wheel position: This forms a linear state constraint such that z s -z us ∈ [z min def , z max def ], where z min def , z max def are the mini- mum/maximum deection position between the chassis and the wheel. Max/Min deection between the chassis and wheel velocity: This forms a linear state constraint such that żs -żus ∈ [ żmin def , żmax def ], where żmin def , żmax def are the mini- mum/maximum deection velocity between the chassis and the wheel.

Figure 3 .Figure 3 . 8 :

 338 Figure 3.7: N s vs J norm CL for dierent values of computational scale factor γ

Figure 3 .

 3 Figure 3.9: dSPACE MicroAutoBox II -1401/1511

Figure 3 .Fig 3 .Figure 3 .

 333 Figure 3.10: Computational eciency between the pNMPC and linearization based MPC controller

Figure 3 .

 3 Figure 3.12: Chassis acceleration for the quarter car system

  proposed MPC controller and the comparison controllers were programmed in MAT-LAB/Simulink environment and was implemented on the INOVE test platform. Two road prole tests were conducted to validate the performance of the proposed MPC controller, which are a) Chirp road prole test and b) Bump road prole test.3.9.2.1 Chirp road prole testThe test involved a chirp road prole with amplitude of 2.5 mm and frequency sweep from 5 Hz to 22 Hz (this corresponds to the comfort frequency range for the INOVE test platform).The road prole is shown in Fig.3.14. The PWM-DC control inputs for dierent controllers is shown in Fig.3.15. It is clearly evident that the proposed MPC utilizes the control authority in a judicious manner such that to minimize the vertical chassis acceleration. The RMS values of the chassis acceleration for the test and the percentage gain with respect to nominal passive damping are listed in Table3.4. The RMS values clearly evinces the fact that the proposed MPC method fares better the nominal passive damping and modied skyhook controller.

Figure 3 .

 3 Figure 3.14: Chirp road prole

  Figure 3.15: PWM-DC input for dierent controllers for chirp road prole

Figure 3 .Figure 3 .

 33 Figure 3.16: Bump road prole

2

 2 Half car mathematical model with stochastic road model Here the objective is to incorporate the road model in to the problem formulation. The mathematical model of the vehicle is comprised of a) Dynamical model of vertical motion of half car vehicle extended with road model and b) Kinematics model for longitudinal motion of the vehicle. Let the left and right corners of the vehicle be indexed with i ∈ {l, r} respectively.

  prole is primarily dependent upon two factors which are a) Longitudinal velocity of the vehicle and b) Road roughness coecient, i.e. the road surface [Tyan et al. 2009]. In accordance with the ISO-8608 standard [ISO 1995], the factored power spectral density (PSD) of road surface is dened with

2

 2 requirements for Parallelized pNMPC without road model The objective design can be briey classied as a) comfort and b) ride handling objective [Savaresi et al. 2010]. 4.3. Parallelized pNMPC scheme for control of semi-active suspension system without road model 71 Comfort objective: The goal of the comfort based objective is to minimize the vertical acceleration of the chassis (z s ), obtained from equation (4.1). The comfort objective for a given look ahead period T l is expressed with : The goal of the ride handling objective is to minimize the roll angle (θ) of the vehicle. The ride handling objective for a given look ahead period T l is expressed with Constraint requirements for Parallelized pNMPC without road modelThe constraints for the SA suspension system primarily arise from the physical limitations and secondarily from performance requirements[START_REF] Baumal | Application of genetic algorithms to the design optimization of an active vehicle suspension system[END_REF]]. For the MPC design, the included constraints are 1. ER-SA damper input constraints: (a) Damper force constraint (Physical): The ER-SA damper force is bounded, i.e. u i ∈ [u min,i , u max,i ], ∀i ∈ {l, r}. (b) PWM-DC input constraints: The operating DC for the PWM signal is constrained to φ i ∈ [φ min,i , φ max,i ], ∀i ∈ {l, r}.

2 .

 2 State constraints: (a) Stroke deection constraint (Physical): This forms a linear state constraint i.e. z d,i ∈ [z dmin,i , z dmax,i ], ∀i ∈ {l, r}. (b) Wheel rebound constraint (Performance): This bounds the deection position between the wheel and road. This is to ensure the tyre deection forces are bounded i.e. z us,i -z r,i ∈ [z rebmin,i , z rebmax,i ], ∀i ∈ {l, r}.(c) Unsprung mass displacement constraint (Performance): This bounds the displacement of the unsprung mass. This reinforces the road holding condition for the vehicle i.e. z us,i ∈ [z usmin,i , z usmax,i ], ∀i ∈ {l, r}.

( a )

 a The rst step is the data initialization where the model/constraint parameters and GPU parameters are initialized. b s , g s represents the size of the grid and blocks respectively[START_REF] Sanders | ghe y exmpleX n introdution to generlE purpose q progrmming[END_REF]. n φ l , n φr represents the number of quantized levels of the PWM-DC input signal for the left/right corner of the vehicle respectively.

4. 3 .

 3 Parallelized pNMPC scheme for control of semi-active suspension system without road model 73 (c) The decorators __HOST__, __GLOBAL__, __DEVICE__ denotes the function call made from host (CPU) to host, host to device (GPU) and device to device respectively. The rst function invoked is the MAIN. 2. MAIN function: (a) The lines 2-3 dynamically allocates memory (DMA) in the host and the device for the objective and inputs for every input combination. (b) In line 8, the KERNEL function is launched with appropriate launch parameters g s and b s . (c) The line 9 transfers the computed objective and input data from the device to the host. (d) The lines 6-8 returns the optimal input U * by nding the index of the minimum objective or constraint violation. 3. KERNEL function: (a) The line 11 sets the thread indices for each PWM-DC input combination. From this point onward, each thread parallely computes the solution for each input combination.(b) The line 12 serializes the 2D grid into a single vector for objective and input vector designation.

4 .

 4 GRID2D function: (a) The line 20-19 sets the PWM-DC input values for left and right corner of the vehicle. (b) The function assigns the input values for every thread call, i.e. for all input comline 25 initializes the objective and constraint violation variable to zero and line 26 invokes the GRID2D function for specic thread indices which corresponds to an input combination.(b) The lines 27-35 runs the Euler integration scheme for the system (4.8) until the look ahead period T l with an integration step of h. In case the constraints are violated, the objective is set to a very high value MAX and the constraint violation is quantied with N function, which computes the 2-norm for the constraints. Otherwise, the objective in(4.11) is computed numerically. (c) The line 36 returns the sum of objective and constraint violation for the given thread indices.Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active suspension system Implementation 1 Parallelized pNMPC implementation Data Initialization: Model/Constraint parameters, b s , g s , n φ l , n φr Input: X 0 , D 0 , Γ 0 ; Output: U * 1: function __host__ MAIN(X 0 , D 0 ,

Figure 4 . 1 :

 41 Figure 4.1: CPU vs GPU pNMPC computation time

4. 4 . 3

 43 Road prole simulation test -Ride handlingThe experiment involved a lopsided bump road prole with an amplitude of 4mm with and the duration of simulation was 10s. The road prole is illustrated in Fig.4.2. Three dierent methods were adopted to analyze the performance of the system which are 1) Nominal damping, 2) ACADO-qpOASES NMPC controller and 3) Parallelized pNMPC method. The settings for each method are listed below 1. Nominal damping (Passive): The PWM-DC input for both left and the right corner was set to a constant value 0.225, i.e. the case when the damper control is switched o and

3 .

 3 Parallelized pNMPC controller: The number of discretization points (complexity) was {n φ l , n φr } = {8, 8} for both left and right corners of the vehicle.

Figure 4 . 2 :

 42 Figure 4.2: Road prole for ride handling

Figure 4 . 3 :

 43 Figure 4.3: Ride handling -roll angle θ

Figure 4 . 4 :

 44 Figure 4.4: Computational time distribution

Figure 4

 4 Figure 4.5: PWM-DC input for Parallelized pNMPC method

Figure 4 .

 4 Figure 4.6: PWM-DC input for ACADO-qpOASES NMPC controller

4. 5 . 5 Figure 4 . 7 :

 5547 Figure 4.7: Histogram of computation time on dierent embedded GPU plaftforms

Figure 4 . 8 :

 48 Figure 4.8: Pareto optimal front between comfort and ride handling objective

  4.11 displays the roll Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active suspension systemangle plot for ISO-E road prole.

Figure 4

 4 Figure 4.9: Longitudinal acceleration prole (a x )

Figure 4 .

 4 Figure 4.10: Longitudinal velocity prole (v x )

Figure 4 .

 4 Figure 4.11: Roll angle (θ) plot for dierent controllers for ISO-E road prole

Figure 4

 4 Figure 4.12: Chassis acceleration (z s ) plot for dierent controllers for ISO-C road prole

hForm 2 -

 2 i (p) ≤ 0, ∀i = {1, . . . n c }(5.5) The two forms of constraint scalarization are Form 1 -The scalar function g can be expressed as a sum over all the maximum of inequality violating constraints (if any) i.e. The scalar function g can be expressed as maximum over all the inequality constraints, i.e. based BBO (Uni-variate case)Consider the optimization problem dened in(5.4) for an uni-variate case, then the optimization problem is dened with min p∈[p min , p max ] J(p) s.t. g(p) ≤ 0(5.8)

Figure 5 . 3 :

 53 Figure 5.3: Graphical interpretation of the terms used in SQP BBO method

  , b f , c f dene the coecients of the approximated quadratic function. As there exists three unknown coecients and three sets of equations dened for the functions values

  based BBO (Multi-variate case)

3 Figure 5 . 4 :

 354 Figure 5.4: pNMPC code generation process

Figure 5 . 5 :

 55 Figure 5.5: Parse tree structure

  as a benchmark S/W to study the performance and computation time of the two examples. The former example was implemented with the assumption of a white box model and the latter example was implemented with the assumption of a black box model. The examples were simulated in MATLAB/Simulink on a Intel Core i7, 16GB RAM PC. The pNMPC C++ codes for the respective examples are listed in Appendix A.

Figure 5 . 6 :Figure 5 . 7 :Figure 5 . 8 :

 565758 Figure 5.6: Cart position, cart force and the pole angle of the system for initial condition x(0) = [0, 0, π 2 , 0] (Case 1)

Fig. 5 .

 5 Fig. 5.8 illustrates the cart force, cart position and the pole angle of the system for

Figure 5 . 9 :

 59 Figure 5.9: Parameterization p 1 , p 2 , p 3 , p 4 , p 5 (Case 2)

Figure 5 .-u max 2 ≤ u 2 ≤ u max 2 -

 5222 Figure 5.10: PVTOL Y-position, Z-position, Roll angle (θ)

Figure 5 .Fig. 5 .

 55 Figure 5.11: PVTOL inputs (u 1 , u 2 )

Fig. 5 .Figure 5 .

 55 Fig. 5.13 illustrates the computation time on dSPACE MABXII. The mean and maximum computation time hovers around 586.04 µs and 745.36 µs respectively and the generated codewas doable under RT conditions. Fig.5.15 illustrates the nonlinear frequency response (refer[START_REF] Savaresi | emiEtive suspension ontrol design for vehiles[END_REF]) from road prole to chassis position for minimum, nominal, maximum and pNMPC controller. It is clearly evident that the pNMPC controller attentuates the road disturbance and performs better than other passive damper settings. Fig.5.14 illustrates the PWM-DC signal, damper force and stroke deection respectively. The PWM-DC signal is chosen judiciously to satisfy both minimization of objective as well as constraint satisfaction.

Figure 5 .Figure 5 .Figure 5 .

 555 Figure 5.14: PWM-DC signal, Damper force and stroke deection plots

  max , θ max and ω max represents the maximum bounds over the input, swing angle and swing angular rate respectively. The OCP for the tracking problem is given as min x(.), p(.)(x(t f ) -x d (t f )) T Q f (x(t f ) -x d (t f )) + t f 0 (x -x d ) T Q(x -x d ) + u T Ru) dt subject to(5.40),(5.41), x(0) = [0, 0, 0, 0] u(p) = p(5.42) where t f , Q, Q f R and T s represents the look ahead period, quadratic stage state cost, terminal state cost, input cost and sampling period respectively. The reference tracking vectorx d = [κ, 0, 0, 0],where κ is the reference signal for position tracking. In order to validate, verify and assess the performance, eciency and computation time of the GPU version of the pNMPC controller, the input parameterization was considered to be of length 128 and subsequently compared against the CPU version of the pNMPC controller. The input parameterization vector is p = [p 1 , p 2 , . . . p 128 ], where the control points are placed equidistantly over the prediction horizon and follows a linear prole. Once the OCP is solved for p * , the input u(p * ) is injected into the system over the period T s and this process is repeated in a receding horizon fashion. The pNMPC S/W settings for both the CPU and GPU version were Integrator -4th order explicit Runge Kutta integrator, Step size -0.2s, Prediction horizon -10s, Number of SQP iterations -8, Constraint form -CONST_FORM::FORM_1 and for the rest, default parameters were utilized.

Fig. 5 .Figure 5 .Fig. 5 .Figure 5 .Figure 5 .

 55555 Fig.5.17 illustrates the crane position, swing angle and the angular velocity for pNMPC-GPU and pNMPC-CPU controllers respectively. As from the plot, the dierence in perfor-

5. 9

 9 Future works and conclusionsThe crux of this chapter is to present a derivative free pNMPC code generation S/W and to validate its performance by means of simulation on multiple examples. From the simulation study as well as from HiL simlations conducted, it is certainly evident that the proposed S/W has applications for several engineering systems, where the model exists either as computer codes or as a data driven model. Despite the S/W is suce for several real world applications, there are certain directions for improvement. The future line of work can be stratied into two parts a) Technical challenges and b) Software challenges. A detailed examination of the aforementioned division of work is expounded below

  PDEs, DAEs and Hybrid systems -In the future, the S/W would include support for control of partial dierential equations (PDEs), dierential algebraic equations (DAEs) and hybrid systems. The overarching goal is to widen the scope of the S/W and provide a one-stop shop pNMPC solution.S/W challenges -The current S/W provides a primitive interface to MAT-LAB/Simulink.In the future, it is planned to provide better interfaces to MAT-LAB/Simulink, Python and Julia to benet all the embedded programmers across the board. D i f f E q u a t i o n d6 = a l p ; END_DIFFERENTIAL / / C o n s t r a i n t s BEGIN_CONSTRAINTS BEGIN_REGULAR_CONSTRAINTS S c a l a r C o n s t r a i n t G1 = { a_l <= a <= a_u } ; S c a l a r C o n s t r a i n t G2 = { a l p _ l <= a l p <= alp_u } ; S c a l a r C o n s t r a i n t G3 = {=15* PI / 1 8 0 <= x5 <= 1 5 * PI / 1 8 0 } ; END_REGULAR_CONSTRAINTS END_CONSTRAINTS / / L a g r a n g i a n and Mayer C o s t BEGIN_OBJECTIVES // L a g r a n g i a n C o s t BEGIN_LAGRANGIAN S c a l a r O b j e c t i v e LC1 = t r a n s p o s e ( Xs ) *Q* ( Xs ) ; S c a l a r O b j e c t i v e LC2 = t r a n s p o s e ( Us ) *R* ( Us ) ; END_LAGRANGIAN / / Mayer C o s t BEGIN_MAYER S c a l a r O b j e c t i v e MC1 = t r a n s p o s e ( Xs ) * Qf * ( Xs ) ; END_MAYER END_OBJECTIVES / / I n p u t p a r a m e t e r i z a t i o n ControlParamZ<L i n e a r >{2 , p1 , a_l , a_u } ; ControlParamZ<L i n e a r >{2 , p2 , a l p _ l , alp_u } ; / / PNMPCGEN s i n g l e t o n o b j e c t PNMPCGEN* pNMPC = PNMPCGEN : : g e t S t o n ( ) ;

  ) ; pNMPC =>setConstForm (CONST_FORM : : FORM_1 ) ; pNMPC =>s e t I n t e g r a t o r (INTEGRATOR : : RK45 ) ; / / G e n e r a t e C c o

Table of

 of

		Acronyms and Notations
	MPC	Model Predictive Control
	OCP	Optimal Control Problem
	pNMPC	Parameterized Nonlinear Model Predictive Control
	SS-pNMPC Scenario Stochastic Parameterized Nonlinear Model Predictive Control
	GPU	Graphic Processing Unit
	EMPHYSIS Embedded systems with physical models in the production code software
	HiL	Hardware in the Loop
	GPGPU	General Purpose computing on Graphic Processing Unit
	FMI	Functional Mockup Interface
	eFMI	Embedded Functional Mockup Interface
	COG	Centre of Gravity
	DAQ	Data Acquistion system
	PCVC	Probabilistic Constraint Violation Certicate
	IP	Intellectual Property
	BBO	Black Box Optimization
	3AC	Three Address Code
	LTI	Linear Time Invariant
	R	Real values set
	C	Complex values set

  1.4. Performance and objective requirements for vehicle suspension control 15

			l l			l r
				m s , I x	θ		z s
						C.O.G
	k s,l	m us,l	u l	z us,r	k s,r	m us,r	u r z us,l
	k t,l	z r,l		k t,r	z r,r

Figure 1.4: Half car roll oriented vehicle model with

  QP ⊂ QCQP ⊂ SOCP ⊂ SDP ⊂ CP. In spite of

	Chapter 2. Theoretical background
	could be classied as derivative based methods and derivative free methods. The derivative
	based optimization methods relies upon the zeroth order (function values) and rst order
	(Jacobian)/second order (Hessian) derivative information of the functions involved in the
	NLP problem. The derivative free methods relies only upon the zeroth order information i.e.
	. Most of the convex problems could be broadly classied into the
	following classes: Linear Programming (LP), Quadratic programming (QP), Quadratically
	Constrained Quadratic Programming (QCQP), Second Order Conic Program (SOCP), Semi-
	Denite Programming (SDP) and Conic Program (CP) as illustrated in Fig. 2.1. It is also
	important to note that each of the aforementioned convex problems form a hierarchy chain
	wherein each class is a subset of another under specic assumptions on the structure of the
	objective and constraints set i.e. LP ⊂ the hierarchy in this classication, the methods involved for solving these convex problems
	are compositionally distinct from each other as the individual solvers for each of these classes
	exploit the structure of the functions or constraints set in a dierent way in order to compute
	the solution eciently within a fewer iterations. It is also important to note that LP is a
	subclass of convex NLP where the objective and constraints sets are linear and polytopic in
	nature and the rest follows suit.
	Nonlinear optimization
	Convex problems Non-Convex problems
	LP QP QCQP SOCP SDP CP
	Figure 2.1: Classication of nonlinear optimization problems
	The problems which are not convex in nature are deemed as non-convex problems, therefore
	by denition only a local minimizer could be obtained and there isn't any guarantee or method
	to establish the global optimum of the obtained solution. There exists variety of methods
	for solving both convex and non-convex problems, however in a broader sense the methods

Table 2

 2

	NLP solver	Method	Interfaces	Programming language
	CVXOPT	IPM	Python	Python
	IPOPT	IPM	AMPL, CUTEr, C, C++, f77	C++
	KNITRO	IPM SQP	AIMMS, AMPL, GAMS C, C++, f77, Java, Excel Mathematica, MATLAB, MPL	C++
	LOQO	IPM	AMPL, C, MATLAB	C
	CONOPT	SQP	AIMMS, GAMS	Fortran
	FilterSQP	SQP	AMPL, CUTEr, f77	Fortran77
	LINDO	SQP	C, MATLAB, LINGO	Proprietary S/W
	LRAMBO	SQP	C	C/C++
	NLPQLP	SQP	C,f77,MATLAB	Fortran77
	NPSOL	SQP	AIMMS, AMPL, GAMS MATLAB, C, C++, f77	Fortran77
	SNOPT	SQP	AIMMS, AMPL, GAMS MATLAB, C, C++, f77	Fortran77
	SQPlab	SQP	MATLAB	MATLAB
	fmincon	IPM SQP	MATLAB	MATLAB
	qpOASES	SQP	MATLAB/Simulink, Octave, SCILAB Python	C++
	FORCES PRO	IPM	MATLAB, Python	C/C++

.1: Second and quasi-second order NLP solvers Derivative free methods can be broadly classied into a) direct search methods, b) model-based methods, c) heuristics based approach and d) randomized algorithms. A detailed exposition into derivative free methods are propounded in the following literatures [Nocedal and Wright

Table 2 .

 2 2: Derivative free optimization solvers

	Solver	Constraints	Bounds	Programming language
	ASA	no	required	C
	BOBYQA	no	required	Fortran
	CMA-ES	no	optional	MATLAB
	COBYLA	yes	required	Fortran
	DAKOTA/DIRECT	yes	required	C++
	DAKOTA/EA	yes	required	C++
	DAKOTA/PATTERN	yes	required	C++
	DAKOTA/SOLIS-WETS	yes	required	C++
	DFO	yes	required	Fortran
	fminsearch	no	no	MATLAB
	GLOBAL	no	required	MATLAB
	GLODS	yes	required	MATLAB
	HOPSPACK	yes	optional	C++
	IMFIL	yes	required	MATLAB
	MCS	no	required	MATLAB
	NEWUOA	no	no	Fortran
	NOMAD	yes	optional	MATLAB
	PSWARM	yes	required	MATLAB
	SID-PSM	yes	optional	MATLAB
	SNOBFIT	no	required	MATLAB
	TOMLAB/GLCCLUSTER	yes	required	MATLAB
	TOMLAB/LGO	yes	required	MATLAB
	TOMLAB/MULTIMIN	yes	required	MATLAB
	TOMLAB/OQNLP	yes	required	MATLAB
	PDF-MPC Package	yes	required	MATLAB

arises from the physical, performance or path constraints imposed on the system and the equality constraints arises from ODE/DAE equations, point constraints and mostly are added for merely the convenience of solving the nonlinear OCP problem via NLP solvers. It is important to note that the non-linear OCP (2.5) is an innite dimensional problem. Nonlinear model predictive control (NMPC) implements the non-linear OCP

(2.5)

in practice by means of receding horizon control technique over the prediction horizon T . The optimization problem in

(2.5

  Computational eciency test . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.8.2 Chirp test with comfort objective . . . . . . . . . . . . . . . . . . . . . . . 58 3.9 Experimental implementation of pNMPC controller on INOVE test platform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.9.1 Comparison controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.9.2 Results and Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 60

	Chapter 3
	Experimental implementation of pNMPC scheme for control for semi-active suspension system
	Introduction to GPGPU Computing for Control Part II pNMPC with RT applications for control of semi-active suspension system
	Engineers using CUDA Programming Paradigm

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.1.1 Related works . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.1.

2 Chapter contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 Control oriented ER semi-active damper modeling and parameter identication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.1 Vehicle modeling -Quarter car model . . . . . . . . . . . . . . . . . . . . 44 3.3 Quasi-static nonlinear ER damper model . . . . . . . . . . . . . . . . . 45 3.4 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.1 ER semi-active damper response time estimation . . . . . . . . . . . . . . 46 3.4.2 Design of experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.3 Non-linear least squares (NLS) based data tting . . . . . . . . . . . . . . 49 3.5 pNMPC design requirements for semi-active suspension system . . . 50 3.5.1 Objective requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5.2 Constraint requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.6 Parameterized NMPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.7 RT HiL implementation of pNMPC controller and Linearization based MPC controller on dSPACE MABXII . . . . . . . . . . . . . . . 53 3.7.1 Linearization based MPC design . . . . . . . . . . . . . . . . . . . . . . . 53 3.7.2 Simulation analysis for pNMPC method . . . . . . . . . . . . . . . . . . . 54 3.8 Real-time Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.8.1 3.10 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 42 Chapter 3. Experimental implementation of pNMPC scheme for control for semi-active suspension system

Table 3

 3

	.1: Estimated ER semi-active damper parameters
	Parameter	Symbol	Value (SI unit)
	Force parameter	f c	21.38(N )
	Deection position parameter	a 1	178.93(1/m)
	Deection velocity parameter	a 2	23.21(s/m)
	Nominal damping coecient	c 0	71.03(N s/m)

Table 3 .

 3 3.7. RT HiL implementation of pNMPC controller and Linearization based MPC controller on dSPACE MABXII 53 2: Model parameters for INOVE quarter car platform and proposed MPC design

	Parameter	Symbol	Value (SI unit)
	Chassis quarter car mass	m s	2.27(kg)
	Unsprung mass	m us	0.25(kg)
	Suspension stiness	k s	1396(N/m)
	Tyre stiness	k t	12270(N/m)
	Max/Min damper force	u, u	±21(N )
	Max/Min deection position	z max def , z min def	±0.005(m)
	Min PWM duty cycle	φ min	0.1
	Max PWM duty cycle	φ max	0.35
	Look ahead period	T l	0.23(s)
	Estimated sampling period for the INOVE platform	T sest	0.023(s)
	Natural sampling period	T s	0.005(s)

3.7 RT HiL implementation of pNMPC controller and Linearization based MPC controller on dSPACE MABXII

Table 3 .

 3 3: RMS values for comfort objective for chirp road prole

	Objective	Linearized MPC	Parameterized NMPC
	Comfort (m/s 2 )	2.7701	2.2643
	3.8.2 Chirp test with comfort objective

Table 3 .

 3

	4: RMS values for comfort objective for chirp road prole
	Controller	RMS (m/s 2 )	% Gain
	Nominal passive damping	6.87	0
	Modied Skyhook controller	6.66	3.05
	Proposed pNMPC controller	6.42	6.5
	3.9.2.2 Bump road prole test
	The INOVE test platform was excited with bump road prole, shown in Fig. 3.16 with peak
	amplitude of 7 mm and duration of 10 s. The recorded chassis acceleration is shown in Fig.
	3.17. From the chassis acceleration plot, it is evident that the proposed MPC method mitigates
	the peak chassis acceleration at bump points. The PWM-DC control input is shown in Fig
	3.18.

5 :

 5 DevicetoHostCpy(h obj , d obj ); DevicetoHostCpy(h inp , d inp ) : function __global__ KERNEL {bs,gs} (d obj , d inp , X 0 , D 0 , Γ 0 ) 11: i ← bIdx.x*bDim.x+tIdx.x; j ← bIdx.y*bDim.y+tIdx.y

	6:	k opt ← indexmin(h obj )
	7:	U * ← h inp (k opt )
	8:	return U *
	9: end function
	12:	r ← n φ l i + j
	13:	if i ≥ n φ l ∨ j ≥ n φr then
	14:	return
	15:	end if
	16:

10

Table 4 .

 4 2 and 4.3 and from the tables, it is evi-

	Chapter 4. GPU based parallelized pNMPC scheme for control of semi-active
	suspension system

Table 4 .

 4

	2: ACADO-qpOASES NMPC controller
	N s FA Mean CT (ms) Max CT (ms) NCLO
	5	0.70	1.6
	10	0.95	2.2
	15	1.3	2.9	0.6484
	20	1.5	3.0	0.6412
	25	1.9	3.9	0.6317
	dent that the proposed parallelized pNMPC method performs better than ACADO-qpOASES
	NMPC controller in all criteria.

Table 4 .

 4 The corresponding best setting was utilized, Integrator -4th order Runge Kutta integrator, QP solver -qpOASES, Hessian approximation -Gauss-Newton, Discretization -Multiple shooting, Discretization intervals -5 and for the rest, default parameters were utilized.

		3: Parallelized pNMPC method
	{n φ l , n φr } FA Mean CT (ms) Max CT (ms) NCLO
	{2, 2}	0.35	0.62	0.4679
	{4, 4}	0.35	0.60	0.4640
	{8, 8}	0.35	0.61	0.4646
	{16, 16}	0.36	0.55	0.4588
	{32, 32}	0.41	0.67	0.4568
	the system is passive.
	2. ACADO-qpOASES NMPC controller:

  Damper force constraint (Physical): The ER-SA damper force is bounded, i.e.u i ∈ [u min,i , u max,i ], ∀i ∈ {l, r}.(b) PWM-DC input constraints: The operating DC for the PWM signal is constrained to φ i ∈ [φ min,i , φ max,i ], ∀i ∈ {l, r}. d,i ∈ [z dmin,i , z dmax,i ], ∀i ∈ {l, r}.

	4.5. Scenario-stochastic pNMPC scheme for control of semi-active suspension
	system	81
	(a) 2. State constraints:
	(a) Stroke deection constraint (Physical): This forms a linear state constraint i.e.
	1. ER-SA damper input constraints:

z

  GPU based parallelized pNMPC scheme for control of semi-active suspension system obtained as return arguments. Finally, the expected objective function EObj i and probabilistic constraint violation certicate PCVC i are numerically obtained for the i th input in U 1:ng .

	Chapter 4.
	84
	:	end for
	29:	Obj ← Obj + TObj
	30: 31: 32: end function end for return { Obj Nτ , CV Nτ }

Table 4 .

 4

	4: NVIDIA boards H/W conguration and results
	Board	#Cores	CC	Mean CT (ms)	Max CT (ms)
	Jetson Nano	128	5.3	20.81	29.83
	Jetson TX1	256	5.3	22.37	30.64
	Jetson TX2	256	6.2	7.62	7.74
	Jetson Xavier	512	7.2	6.04	6.59

Table 4 .

 4

		5: RMS value of chassis acceleration
	ISO	SS-pNMPC(ms -2 )	u n (ms -2 ) u(ms -2 ) u(ms -2 )
	A	0.012	0.014	0.011	0.017
	B	0.023	0.029	0.021	0.034
	C	0.049	0.058	0.042	0.076
	D	0.093	0.118	0.087	0.138
	E	0.194	0.234	0.178	0.277

  which is based on indirect methods and the toolbox generates C code of a RT algorithm to update the initial co-state of the OCP problem. A newer version of the toolbox AutoGenU [Ohtsuka 2015] provides extension with continuation/GMRES method and a RT optimization (RTO) algorithm with interface to Maple. Ever since AutoGen, there has been signicant development in building code generation based MPC toolbox for RT applications.

  = [p 1 (t), p 2 (t), p 3 (t)] and

	100										Chapter 5. pNMPC code generation tool
		1.4				Parameterized sinuoid input
		1.2							Parameterized sinusoid input
				1
	Sinusoid input	0 0.2 0.4 0.6 0.8
		-0.2
		-0.4
		-0.6
				0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
											Time(s)
					Figure 5.2: Parameterized sinusoidal control input (u(p(t), t))
				1.5				Amplitude, frequency and phase of parameterized sinuoid input
		Amplitude	-0.5 0 0.5 1	(t 0 ,p 1,0 )	(t 1 ,p 1,1 )	(t 2 ,p 1,2 )		(t 3 ,p 1,3 )	Parameterized amplitude (t 4 ,p 1,4 )
				-1
					0	0.2		0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
											Time(s)
		Frequency(Hz)	-0.2 0 0.2 0.4 0.6	(t 0 ,p 2,0 ) (t 1 ,p 2,1 )		(t 2 ,p 2,2 )	(t 3 ,p 2,3 )	Parameterized frequency (t 4 ,p 2,4 ) (t 5 ,p 2,5 )
			-0.4
					0	0.2		0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
											Time(s)
				1.5
		Phase	0.5 1	(t 0 ,p 3,0 )		(t 1 ,p 3,1 ) (t 2 ,p 3,2 )	(t 3 ,p 3,3 )	(t 4 ,p 3,4 )	(t 5 ,p 3,5 )	Parameterized phase (t 6 ,p 3,6 ) (t 7 ,p 3,7 )
				0
					0	0.2		0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
											Time(s)
	Figure 5.1: Parameterized amplitude (p 1 (t)), frequency (p 2 (t)) and phase (p 3 (t)) for sinusoidal
	control input

Table 5 .

 5

		1: Notation and meaning
	Notation	Meaning

Chapter 5. pNMPC code generation tool input and the last two arguments represents the minimum and maximum bounds over the input p1 respectively.

ControlParamZ<Linear> -This class parameterizes the input in a linear fashion over the prediction horizon. The other specications follow suit as the piecewise parameterization. Usage: ControlParamZ<Linear> {2,p2,-1,1};

Real/Symbolic classes

The pNMPC software supports several functions and operation between real numbers and as well as symbolic objects. The list of supported functions are {sin, cos, tan, sinh, cosh, tanh, exp, log, abs, asin, acos, atan, asinh, acosh, atanh, minimum, maximum, sign}. The list of supported numerical operations are {+, -, * , \, ∧, ≥, ≤, &&, ||}, where the last two are logical AND and OR operations, which is used to couple multiple constraints to one (Example -blocking constraints in state space region). The scalar real numbers are declared with Real keyword and real valued matrices are declared with MATReal keyword. Symbolic matrices are declared with MATHyperStates keyword. Usage:

MATHyperStates A(2,2); MATReal B (2,2);

There are several in-build matrix operations such as matrix-vector operations, matrix-matrix operations etc. packaged along with the S/W. // qet ingleton instne of ode genertion lss PNMPCGEN* pNMPC = PNMPCGEN::getSton(); pNMPC->setInitialTime(0); pNMPC->setFinalTime(2); pNMPC->setStepSize(0.1); pNMPC->getSolver()->setNiter(4); pNMPC->setConstForm(cF); pNMPC->setIntegrator(iODE); // qenerte g odes pNMPC->genCCode();

Code generation classes

pNMPC code generation module

The pNMPC code generation module is not dissimilar to any compiler design paradigm [START_REF] Aho | Compilers, principles, techniques[END_REF] at the same time the steps involved are not as extensive as for compilation process of any programming language. The motivation for adopting a scheme as such are two folds which are:

1. The ingredients of the OCP problem (objectives, constraints, dynamics etc.) fed by the user are broken down into fundamental elements and then modeled in an appropriate way to suit the optimization module. Case in point, the inequality constraints ought to be aligned in a non-positive formulation as described in equation (5.4).

2. By breaking down the OCP into it's fundamental elements, code optimization can be performed eciently which in turn benets in reducing the memory footprint (space complexity) and burning less computer clock cycles (time complexity). This feature has high practical importance, especially for low-end embedded devices. #d e f i n e PI 3 . 1 4 1 6 u s i n g n a m e s p a c e pNMPC; i n t main ( ) { / / S t a t e s S t a t e s p , t h e t a , pd , t h e t a d ; / / I n p u t s I n p u t s p1 , p2 , p3 , p4 , p5 ; / / C o n s t a n t p a r a m e t e r s R e a l m1 = 1 , m2 = 0 . 1 , g = 9 . 8 1 , l = 0 . 5 ; / / Bounds R e a l pmax = 2 , umax = 1 0 , t h e t a m a x = 2* PI ; / / OCP d a t a MATHyperStates Xs ( 4 , 1 ) ; Xs = {p , t h e t a , pd , t h e t a d } ;

MATReal Q = d i a g ( { 5 , 1 0 , 1 , 1 } ) ; MATReal Qf = d i a g ( { 1 0 , 2 0 , 1 , 1 } ) ; R e a l R = 0 . 1 ; / / I n p u t p a r a m e t e r i z a t i o n BEGIN_PARAMETERIZATION_MAP P a r a m e t e r i z a t i o n M a p u = p1 *p + p2 * t h e t a + p3 *pd + p4 * t h e t a d + p5 ; END_PARAMETERIZATION_MAP / / D i f f e r e n t i a l e q u a t i o n s / / E x t e r n a l v a r i a b l e s E x t e r n a l p s i _ y = " ModelY " ; E x t e r n a l p s i _ z = " ModelZ " ; / / C o n s t a n t p a r a m e t e r s R e a l a_l = 0 , a_u = 2 , a l p _ l = =2, alp_u = 2 ; / / D i f f e r e n t i a l e q u a t i o n s pNMPC =>s e t I n i t i a l T i m e ( 0 ) ; pNMPC =>s e t F i n a l T i m e ( 2 ) ; pNMPC =>s e t S t e p S i z e ( 0 . 1 ) ; pNMPC =>g e t S o l v e r ()=> s e t N i t e r ( 4 / / D i f f e r e n t i a l e q u a t i o n s BEGIN_DIFFERENTIAL D i f f E q u a t i o n d1 = x2 ; D i f f E q u a t i o n d2 = p s i _ y ; D i f f E q u a t i o n d3 = x4 ; D i f f E q u a t i o n d4 = p s i _ z ; D i f f E q u a t i o n d5 = x6 ;