Model-based design (MBD) is gaining popularity in various applications, such as medical, industrial, aerospace, and automotive. This methodology is also sometimes used to develop control systems for guided missiles. However, in the existing literature, finding a detailed and comprehensive description of the process is challenging. Such data are often classified due to security reasons. Moreover, the development of control systems for high-speed objects is very costly and restricted by military institutions. A lot of the existing systems are based on legacy solutions that were developed several years ago. Also, testing the developed system requires specialized test ranges because failure during real flight might have catas-trophic consequences and be very dangerous.
When a malfunction occurs, the object can hit the unintended target accidentally. For this reason, using models plays a critical role in preparing flight tests. The onboard software must be developed cost-efficiently and ensure a final high-quality product. There exists a significant research gap in publications relating to the topic of automatic code generation (ACG) in the field of guided HSUAVs. The overall process workflow is generally well-known, but from existing available literature, itis hard to conclude about the details (for example, settings of the code generator or verification method).
Alpaslan [3] analyzed the challenges with the development of weapon software. The autopilot might be considered as a safety-critical system operating in real-time. The development of such a system is often more complicated than that of commercial products. Also, information security must be considered. Today, a lot of weapon functionalities are included in the software. As a result, requirements can be pretty similar to those in the space industry.
The traditional process of software development was manual coding [7]. However, such an approach is time-consuming. Up to this time, many papers on code generation have appeared in the literature. A review of recent studies on automated code generation was presented by Kshirsagar et al. [39]. ACG was used for the aerospace projects. Maroli et al. [46] adopted code generation for the motor controller for NASA’s X-57 Maxwell aircraft. They also explained the role of several analysis tools available in MATLAB to ensure that the created code is bug-free. This technique was applied to develop the software for the fault protection system in the Deep Space 1 mission [47]. MBD and Simulink were successfully used to create the Guidance, Navigation, and Control (GNC) software for the Orion project, and NASA published a detailed description of the methodology in [31,64]. Fraticelli [23], presented the example of software development for an Attitude Determination Control System (ADCS) of a satellite and investigated the numerical equivalency between the model and C code output. Carpenter [10] also reported using ACG to obtain the ready-to-use software for nanosatellite in the context of ADCS. Yahyaabadi et al. [71] studied the feasibility of using auto-code generators in future space missions.
Arm et al. [5] discussed the code generation problem for safety-critical applications and presented an example of a complete workflow for the ARM Cortex R hardware target. Erkinnen [17] also addressed the topic of software development by using ACG for such systems. Several years ago, Schwarz [59] mentioned that code generation is not often applied for safety-critical systems because such applications should obey strict requirements. However, up to now, auto-coding using automatic code generators is still not a widely accepted method of software creation in safety-critical applications. The main reason behind this is the lack of certification and qualification of the generators [15,22,62]. Moreover, many of the existing aerospace and defense systems were developed a couple of years ago, and there is a problem with adequately integrating the legacy codes with modern software development tools.
Autogenerated code from Simulink can have bugs that might have fatal consequences [34]. The resulting code must work in the same manner as the original model of the system. The problem of numerical equivalency between the model and autogenerated code was addressed by several researchers [6]. Lambersky [40] presented the development of the control algorithm for a legged robot and compared the results from the simulation and target platform.
Many examples of using ACG were reported for relatively simple systems that do not have to meet strict requirements [41,48,52,53,60]. For instance, Hyl and Wagnerová [28] presented the case of the controller development process for the climate unit laboratory stand. Otava [55] showed motor control development using Simulink Embedded Coder and verified the results experimentally. Fakih and Warsitz [21] reported the code generation process for two simple user case studies and compared the results of MIL and SIL simulations. Krizan et al. [38] presented a detailed description of a model-based design for control of a direct current brushless motor. Also, Wu [69] reported using MATLAB built-in tools from the Coder family for controller design for an electric motor using TI FI28379D hardware. Lixin and Liwen [11] presented the case of code generation for an autonomous underwater robot.
Existing works are related to using MATLAB tools in the MBD process. It must be mentioned that a set of other code generators exists. SCADE environment with a qualifiable KCG Code Generator is used for the development of avionics software [67]. In the auto-motive industry, the TargetLink generator reached significant popularity. Ajwad [2] compared the performance of Embedded Coder and TargetLink. Also, Scilab-Scicos and GeneAuto were proposed as open-source solutions [37,58,66]. Jacobitz and Xiaobo proposed [32,33] LoRra generator that can be used to translate Scilab/Xcos models into usable C code. Also, a group of research papers proposes customized, much less popular generators other than widely used (both commercial and open-source) solutions.
For example, Yu et al. [72] designed a custom generator named Mercury and compared the per-formance with Simulink Embedded Coder. Bourbouh et al. [8] proposed a framework named CoCoSim that can be used to analyze and verify the code generated from Simulink models. Some recent studies concentrate on generating parallelized code from Simulink block diagrams [70].
Several trials by many individual researchers and various organizations were performed to standardize the flight-ready production code generation process from Simulink models. For example, Erkkinen [16] described several available guidelines and best practices for automatic code generation for flight applications to make it optimized and efficient. Some years later, in [18], he also presented the use of a qualified verification tool for the model-based designs for DO-178B applications. Fraticelli [24] described a set of valuable practices in generating C code from the Simulink model. MathWorks Automotive Advisory Board (MAAB) proposed a set of proper practices and recommendations. Motor Industry Software Reliability Association (MISRA) initially created a set of C standards only for automotive applications, but now they are widely used in other branches of industry. NASA has its standards and regulations (for example, NPR 7150.2, NASA-STD-8739.8, and NASA-GB-8719.13). One of the most important contributions and progress in that field in the context of aerospace applications is reports prepared by the Space AVionics Open Interface aRchitecture (SAVOIR) working group that can be used as guidelines in C code generation from Simulink models [19,20].
In the field of guided UAVs, the literature related to MBD is relatively scarce. Among the detailed publications, one might mention the work of Putro and Septanto [57], who developed a real-time test environment for evaluating the missile’s performance. Craft [13,14] presented the case of using automatic code generation for the Long Range Land Attack Projectile (LRLAP). Abdelaty [1] gave a comprehensive description of the autopilot development for the anti-tank round. Tancredi [65] described the overall design methodology, and reported using code generation by MBDA company but without describing a particular case. Holliday [27] reported using MBD by Thales company and presented several user cases, however, without sufficient detail. Several authors briefly described Hardware-in-the-Loop simulators for guided high-speed objects [4,36,61,76], but the results are difficult to reproduce. Waxenegger [68] presented a Hardware-in-the-Loop testing methodology for a rocket propulsion control system.
Automatic code generation from the system model offers several advantages compared to the software’s traditional manual coding for the control system. The overall cost and time of the prototyping might be reduced significantly [74]. As a result, auto coding can reduce the overall effort compared to manual coding [51].
The research might concentrate on testing various variants of the algorithm instead of implementing low-level language details [45,75]. The resulting code has a clear and repeatable structure. This technique also decreases the need for manual debugging, which is often time-consuming. Such methodology reduces the probability of introducing human-made mistakes during manual coding (for example, in signs). Fur-thermore, the consecutive versions of the software might be easily managed. The verification and documentation process could also be automated. Auto-matically generated code can be as fast or even faster than handwritten code [9,12,44]. Using MBD, the number of expensive and time-consuming flight trials that must be performed to validate the control system can be reduced significantly. SIMULINK diagrams are quite easy for managers and system engineers to understand, so team collaboration is simplified. On the other hand, this methodology requires the staff to have multidisciplinary competencies (flight dynamics, modeling, electronics, etc.).
The motivation for the presented study was strictly practical. There was a need to develop a control system prototype for the gasodynamically controlled HSUAV in a short time. The most common actuation method is using movable aerodynamic fins or thrust vectoring. Using multiple solid propellant lateral motors is much less prevalent.
Up to now, only a few existing solutions have used this type of actuation. Such pulse thrusters (as only actuators) were used on the M47 Dragon anti-tank and STRIX mortar rounds. Ukrainian Vilkha adopts 90 small pulse thrusters in the front of the fuselage, but also embraces the aerodynamic movable fins used in the last phase of flight. Moreover, 180 lateral motors are used on PAC-3 MSE in the initial and terminal phases of flight (this platform also uses aerodynamic control).
The main contribution of this study is the detailed description of automatic C code generation from the Simulink model for the HSUAV autopilot. The originality of this research lies in the fact that the object is actuated only by solid propellant lateral motors, which is not a standard solution. The control systems for such a configuration are not described in sufficient detail. For a tactical-grade system, the reliability of the operation is a critical factor. The solution must be carefully tested for various possible scenarios that might appear in the actual firings. The resulting code cannot include unintended functionality affecting the system’s operation. This research extends the results presented in [30, 42]. The proposed approach might be helpful for other researchers in the field of external ballistics. The article also includes a description of a set of good code-generating practices.
The structure of the paper is as follows. At first, a description of the flying platform is shown. Then, a flight simulation model of the HSUAV developed in Simulink was presented. Next, the control algorithm was described. Later, the code generation process was explained.
The representative results were shown. The paper ends with a discussion of the obtained results and a summary of the main findings. Finally, some further possible research directions are suggested.
The diameter of the fuselage is 0.122 m and the total length of the object is 2.1 m. At launch, the mass is around 26 kg, and moments of inertia are 0.13 kgm2 (longitudinal) and 10.2 kgm2 (transversal). The mass of the propellant is 6 kg. The operation time of the main motor is a little more than 2.8 s, and the available thrust has a maximum value of about 8600 N. HSUAV is aerodynamically stabilized by four trapezoidal fins (there is no wrap-around functionality). These stabilizers are canted, and as a result, they create a rolling moment (the platform intentionally spins slowly during the flight). The platform is equipped with 32 solid propellant lateral thrusters located in front of the center of mass (no movable aerodynamic fins were used). The thrusters are arranged into four layers with eight motors in each. Each of these thrusters might be used only once. The operation time of the thruster is less than 0.05 s, and the mean thrust is more than 700 N. This kind of actuation is challenging because the object possesses very low control authority. That means the single pulse thruster can translate the impact point location by only several meters. The maximum range is approximately 11 km. In the considered version, the solution accounts for only stationary targets.
The GNC system comprises a navigation unit and an actuation unit. These subsystems are integrated into a control unit using the mechanical structure.
The main components of the navigation unit are a set of sensors, a Primary Control Computer (PCC), and three antennas. The main function of that system is to use data collected by the sensors and convert them into information about the object’s flight state. The navigation data are obtained using an off-the-shelf available Inertial Measurement Unit (IMU), a custom-developed GALILEO receiver, a pressure sensor, and infrared sensors (these sensors are equally spaced around the circumference of the object body). The PCC uses an ARM Cortex A7 processor. The data obtained from all sensors is fused using a Kalman filter. One of the main difficulties of the adopted solution lies in the fact that the roll angle must be estimated very precisely (with accuracy < 1°). The low-cost micro-electromechanical systems (MEMS) gyroscopes that are available on the market have a low measurement range (typically up to 2000°/s) and are not entirely suitable for the developed system (sensor saturation could occur). Due to this, satellite signals needed to be received by the antennas mounted on the sides of the fuselage to estimate the roll angle. Moreover, the infrared sensors are used to improve the calculation accuracy. Additionally, the information about the altitude is corrected by the readings from the pressure sensor.
The actuation unit comprises a Lateral Motors Control Computer (LMCC), a set of power amplifiers, and lateral motors. The primary function of this unit is to use the data about the current platform state (obtained from the navigation module), calculate the steering commands, and activate the appropriate lateral thrusters. The LMCC is based on the ARM Cortex M4 architecture. It calculates when each motor should be fired to steer the HSUAV on target. The power amplifiers activate the igniters attached to the lateral motors. They convert the output data from the LMCC into a set of electrical signals for the 32 ignition channels. Each igniter (sometimes called a motor starter) is responsible for starting the combustion of a single small solid propellant in the lateral motor. This unit also includes an Electrostatic Discharge (ESD) protection circuit to prevent each igniter from unintentional release.
Additionally, the HSUAV is equipped with a radio telemetry downlink that transmits the data (linear accelerations, angular rates, position, etc.) to the ground control station (GCS) with a frequency of up to 500 Hz. The object state can be observed in realtime (in the form of graphs and visualization) on the GCS. During the launch phase, the electronic units are subjected to high acceleration. For that reason, triple-redundant communication channels were used to increase the system’s reliability and minimize the probability of failure in flight.
The onboard battery powers all electronic subsystems described above.
In the MBD approach, one of the first steps is to create a realistic plant model. To make the real prototype of the hardware, the detailed flight simulation model of the real object was developed and implemented in the MATLAB/Simulink environment. The baseline mathematical formulation can be found in “Appendix A”. The detailed description of the mathematical model and its implementation can be found in [26, 29, 30, 42, 43, 56].
The vehicle subsystems (main motor, lateral thrusters, onboard sensors, etc.) were modeled based on the data of real hardware. The dynamics model parameters were obtained mainly by the experiments in laboratory conditions and firings of the existing prototypes (unguided ones). Additionally, they were confirmed using CAD models of the structure. HSUAV moments of inertia were measured using bifiliar and quadfilar torsional pendulums. The thrust curves of the main motor and lateral thrusters were acquired by measurements on the test stand at various temperatures. Aerodynamic parameters were collected by using specialized semi-empirical and Computational Fluid Dynamic (CFD) codes. Later, after initial flight trials, the aerodynamic database was corrected by the appropriate form factors. The sensors’ data were found using datasheets from the manufacturers and confirmed by laboratory experiments.
Also, the launcher dynamics were included in the model to make the simulation of the initial phase of flight more realistic and include several important phenomena (for example, the tip-off effect). Environmental conditions (e.g., air temperature and humidity) are taken into account to predict the atmospheric properties.
The key element of the simulation is the model of the control system hardware and the software. The inputs to the control algorithm are current position coordinates, Tait-Bryan angles (yaw, pitch, roll), angular velocities, linear velocities, and linear accelerations. The output is the firing command of each thruster. The target position coordinates (latitude, longitude, and altitude) must be known and programmed before the launch.
The scheme of the overall control algorithm is presented in Figure 1.
Top-level Simulink scheme of the control algorithm
The quality of the resulting C code can be affected by a lot of factors and settings. To appropriately choose the compiler parameters, the model was created using SAVOIR guidelines [19, 20]. The equations of motion were integrated numerically using a fixed-step Bogacki-Shampine (third-order) solver with a step size of 0.0001 s. Only fixed-step solvers are currently supported in MATLAB during the code generation process. The model cannot contain any blocks that are not supported by the code generator. Before performing ACG, it is crucial to ensure a well-styled Simulink model. All of the inputs and outputs were named using strict predefined conventions. The names and units were displayed on the signal wires to eliminate the probability of human error. Signals were grouped using buses, which ensured readability. Signal line crossing was eliminated. The scheme also includes a set of diagnostic outputs that can be used in later stages to analyze in detail the results on the actual hardware. Equality comparison using floating point numbers was avoided because that can lead to entirely unpredictable results.
A small team of experienced programmers developed the model, so it was necessary to perform version control. Each subsequent version of the code has a unique name with data and time of last saving included. Simulink Design Verifier was also used to detect potential problems with the block diagram (for example, dividing by zero or dead logic).
The requirements and acceptance criteria for individual subsystems were formulated at the initial phase of the control system development. These requirements were also created for the software. The code should be efficient and able to operate in real-time on the target hardware. The control algorithm was carefully optimized to increase the execution speed and save the available computational resources. The autopilot code should be integrated with other software parts created by different cooperators. It is obvious that the software should include appropriate comments and be easy to understand.
The GNC module can be considered as a complicated multi-rate digital system. The data about the individual HSUAV flight parameters can be delivered from the navigation system module with various (often relatively low) frequencies. For example, the position using a GALILEO receiver can be obtained with a frequency of no more than several hertz. On the other hand, the raw data from IMU outputs (angular rates and linear accelerations) could be sampled with a frequency of up to 200 Hz. These data are fused using a Kalman filter [73] to obtain a more accurate estimation of the flight parameters (mainly velocity, position, and Euler angles), but the output frequency is no more than 40 Hz. The thrusters can be fired only in certain phases of flight when the roll rate of the object is approximately 2 revolutions per second (after apogee). If the roll angle is estimated at 40 Hz, then the motors can be fired with the angular resolution 18° which is too low. That means the control loop should operate at a much higher frequency than the data can be delivered from sensors. It was estimated that the control algorithm should be recalculated with a frequency of 5000 Hz. That way, the firing conditions will be checked with an angular resolution of 0.144°.
For this reason, an upsampling algorithm should be applied to obtain the navigation data with the required frequency. On the other hand, such an algorithm cannot introduce time delays because then the object might be steered in the inappropriate direction.
A custom module (whose primary function is to extrapolate the data gathered from the various onboard sensors) was developed to overcome that difficulty. This module consists of two parts: extrapolation of roll angle and position, which were determined to be the most influential on the algorithm’s accuracy. On the output of both extrapolation modules, the navigational information (roll angle and three-dimensional position in space) is delivered with the desired frequency of 5000 Hz. All other signals are obtained directly from the navigation module without the extrapolation procedure (their frequency is sufficient).
The roll angle extrapolation block scheme is presented in Figure 2.
Roll angle extrapolation block scheme (signal frequencies denoted in blue font)
The main idea behind this method is that the angular rate (obtained directly from gyroscopes) can be delivered with a higher frequency when compared to information about the actual roll angle. It was assumed that between two consecutive samples of the roll angle, the information on the intermediate points can be obtained using linear interpolation. The forward Euler method was used for numerical integration of the roll rate to obtain the roll angle Φ:
The same idea was used to extrapolate the position coordinates:
A modified Proportional Navigation Guidance (mPNG) algorithm was developed to steer the object accurately towards the target. The target location
The components of the distance between the target and the vehicle (along each axis of the NED frame) are:
In a similar way, the components of the relative velocities are:
The LOS angles in each plane are:
These angles must be differentiated with respect to time to obtain LOS rates.
Next, the acceleration obtained from equation 5 must be converted to the body-fixed frame and compen-sated by the gravity:
At the outputs, this algorithm delivers information about the commanded linear accelerations in the body-fixed frame. The three components of the accel-eration vector are used as inputs to the Firing Logic (FL) submodule (please see the next section of the paper). To hit the target successfully, the vehicle must realize this commanded acceleration. The commanded direction of flight in the body-fixed frame is calculated as:
The magnitude of the commanded acceleration is:
The details about the mPNG algorithm are presented in [30,56]. The proposed algorithm can achieve a significant impact point dispersion reduction, which was proved by various studies conducted in e.g. [25,26,29,35,63], where very similar types of guidance method were utilized.
Additionally, it must be remembered that the navigation system estimates the velocity
The discrete nature of actuation makes the guidance process difficult because there is no possibility of continuously tracking the commanded lateral accelerations. Since only 32 thrusters are available, their use during the flight must be carefully planned. If the motors are fired too early, then it will be impossible to influence the trajectory in the terminal phase of flight. Moreover, if the time between two pulses is too short, the platform can achieve high angles of attack and sideslip, and as a result, disintegration can occur due to high structural loads. On the other hand, if the thrusters are fired at too large time intervals, their usage will be ineffective, and the HSUAV will not achieve its destination.
The firing logic is based on the conjunction of four conditions that could be expressed mathematically as (for details, please see [30]):
The individual thruster is fired when the output signal from the algorithm changes from 0 (low state) to 1 (high state), which happens after all four conditions are true simultaneously.
The state machine (Figure 3) was also used to realize a step-by-step control process to prevent too-early lateral motor firings for safety reasons. The actuation unit with the control motors might be activated only after detecting the launch event (longitudinal acceleration of more than 50 m/s2 over 0.4 s).
Simplified sequence of events before activation of the thrusters
The lateral thrusters are fired in the predefined sequence prepared as a firing table (block “Motor-FireMatrix_simple” in Figure 1. The system user might program this sequence before the launch. An example of such a table is presented in Table 1, where the first column is the number of the pulse, the second column is the ID of the thruster, and the last column is the angular position of the thruster on the fuselage’s circumference.
Example of the firing table.
| Pulse no. | Thruster ID | Angular position |
|---|---|---|
| 1 | 2 | 0° |
| 2 | 5 | 25° |
| … | … | … |
| 32 | 28 | 285° |
The software should be developed using a strict strategy. Several management methodologies exist that can be applied, but some are unsuitable in this context. For example, the waterfall methodology was not feasible because the initial requirements were not fully formulated at the beginning of the process. The software development process diagram is presented in Figure 4.
Software development process
The onboard hardware is quite limited in computational resources, and it is evident that the control algorithm prototype developed in MATLAB must be implemented in the appropriate low-level programming language. The Simulink Coder (earlier Real Time Workshop) and Embedded Coder were used for automatic code generation. However, this MATLAB extension was insufficient to automate the whole process entirely. For this reason, a custom code was developed to generate, compile, run, and test the algorithm. Microsoft Visual Studio (MS VS) 2019 was used to compile the code.
To generate the code appropriately, the settings must be carefully chosen, as the default code generation settings are not optimal [50]. ISO/IEC 9899:1999 C language standard was used. The code was generated with various settings, and the results were assessed for each combination of parameters. For the generated code, a Visual C\C++ file for Embedded Coder was chosen. Also, the influence of code optimization settings on the accuracy of the generated code was investigated. After a set of trials, it was decided to set the optimization level as a minimum. The algorithm parameters (e.g., target location can be set as inlined or tunable. It was decided that they should be inlined in the code to reduce memory usage and increase code efficiency. In that way, there is no need to preallocate the memory for the parameters.
The embedded Coder Support Package for ARM Cortex-M Processors was used to ensure that the code would be effective in realizing math operations using the Cortex Microcontroller Software Interface Standard (CMSIS) library. Using this addon, it is possible to replace some functions on the more effective implementations (for example, “sin” trigonometric function is replaced on “arm_sin_f32”).
Single-precision float-type numeric numbers were used to save the available memory. By default, MAT-LAB operates on double-precision numbers, but using them in the microcontroller is not necessary and seriously slows down the computations. For that reason, all input and output ports in the Simulink blocks of the control system model were manually set to “single”. In that way, it was ensured that there would be no unintentional conversion between “single” and “double”. It was detected that this issue is crucial and might lead to catastrophic consequences.
Before being used on real hardware, the algorithm must be tested for various operating conditions. A custom test procedure was developed to verify the control algorithm’s reliability and ensure high code coverage.
The code generation process started with the selection of an appropriate flight scenario. At first, a vast number of Model-in-the-Loop (MIL) simulations (for individual flight scenarios) were evaluated to detect potential problems at the early stage of control system development. At this stage, only a control module model was connected with the HSUAV model. These tests were realized for various launcher azimuth and elevation angles, various thrust curves of the lateral thrusters, etc. Next, the structure of the Simulink control system model was iteratively optimized to achieve the best possible configuration and decrease execution time. The model was carefully checked to prevent inefficient parts of the code (for example, algebraic loops). Not all functions are supported for code generation. In that case, a custom equivalent block structure replaced such an unsup-ported function. Simulink Profiler was used to measure the execution time of individual model components. Simulink “Accelerator mode” was set to speed up the execution time. In that way, the resulting time of a single simulation was decreased to only 7 s (for 50 s of the real flight). Then, Monte Carlo simulations were used extensively to evaluate accuracy and precision. The Parallel Computing Toolbox was used to assess the simulations on the workstation (Intel® Core™ i9 and 16 GB RAM). In that way, it was possible to simulate thousands of firings and estimate the measures of dispersion (for example, Circular Error Probable). Also, parametric analysis was conducted to investigate the influence of individual autopilot settings on the overall accuracy and to understand the system behavior. The example of results from this methodology can be found in [26]. This phase ends when the controller model meets the predefined requirements and generates appropriate outputs.
After the MIL stage, detailed documentation was created to describe the low-level functionality of the flight simulation. The documentation of the results was based on the automatic generation of structured reports. After each simulation in MATLAB, two text files were recorded and saved on the computer disc. The first file included in the tabular form time and platform flight parameters: angular rates, position, Euler angles, and linear accelerations. The second file included the array of firing commands. In the later stages of the whole testing process, the abovementioned two files were used as reference data (in the end, the output of the final code should match the reference firing commands).
After completing the MIL phase, the next step was to perform Software-in-the-Loop (SIL) tests. At the SIL stage, the C code from only the autopilot module was generated and used in the control loop instead of the controller system model. The traceability function that is available in Embedded Coder was used extensively to navigate between the Simulink model and C code. These tests were performed on the standard desktop computer (Intel® Core™ i7, 32 GB RAM, MS Windows 11 operating system) without dedicated hardware. At this stage, checking if the autopilot software could be implemented was possible.
Additionally, the code was tested externally to ensure no problems with the generated C software. The above-mentioned custom-developed test procedure was as follows. The C code was run on a Desktop computer in MS VS, and the results (lateral thruster’s firing commands) were stored in a separate text file. Then, both text files (reference data from MATLAB obtained at the MIL stage and results from the standalone C code) were compared. If the results met the acceptance criteria, the C code was ready to be implemented on the real hardware and tested again.
If any anomaly output or unintended behavior was observed, then it was necessary to return to the MIL stage, introduce modifications, confirm that after adjustments, the MIL results were still acceptable, and then again continue the SIL level using the modified software. This process was automatized using MAT-LAB scripts.
Next, the Processor-in-the-Loop (PIL) phase was performed. The developed algorithm in the form of C code was implemented on real hardware with an ARM Cortex M4 processor. The developed control system was tested in laboratory conditions. At first, the (PIL) simulation was evaluated not in real-time. The reference data, including flight parameters, were sent from the text file to the hardware by the interface (these data were used instead of the information from the navigation unit). The control computer calculated steering commands, and the outputs from the system were logged into text files. The main goal of that stage was to check that the results were the same as the reference data. Then, real-time calculations were performed to ensure that the control loop was able to operate with a sufficient frequency (5000 Hz, as mentioned before).
One of the last steps of the testing was the Hardware-in-the-Loop simulation. At this stage, the controller performance can be verified (at least partially) without using the complete real plant. This approach is much safer than testing the control system on the real object in flight tests. The system operated in real-time.
In the last stage, the results from MIL, SIL, PIL, and HIL were compared to ensure that the results were in agreement. If the results of HIL were unsatisfactory, the process came back to MIL, SIL, or PIL levels. Otherwise, if the code met the acceptance criteria, it was decided that the software might be implemented on the real vehicle.
The main goal of the experiments was to check the numerical equivalency (or eventually similarity) of the results between the model in Simulink and the resulting C code that might be directly implemented on the real hardware. The obtained C code must work like the original MATLAB program when tested with the same input data. The numerical discrepancy can be caused by various issues, for example, low-level implementations of the trigonometric functions in the two programming languages [54].
MATLAB R2023 with Update 5 was used to perform numerical simulations. The control algorithm was improved step by step, iteratively. Several iterations of MIL, SIL, and PIL simulations were evaluated to meet the acceptance criteria. The solution was acceptable if the impact point dispersion measured using Circular Error Probable 50% (CEP50%) was smaller than 25 m and the results were repeatable for various flight scenarios. Here, the selected example of a mission scenario is presented.
The following initial conditions were assumed: launcher elevation angle 45.2° and azimuth 46.5°. The algorithm settings were as follows: tt = 0.25 s (minimum time between two pulses), Θg = 5° (pitch angle threshold), tg = 30 s (threshold time), γt = 1.5° (angular tolerance of lateral motor firings), at = 0.5 m/s2 (threshold value of the lateral acceleration). The launcher was located at the origin of the NED coordinate system. It was assumed that the target is stationary, and its coordinates in the NED frame were set to (6430.7, 6222.0, 0.0) m.
At first, MIL tests were evaluated for two cases: unguided flight (the control system was intentionally deactivated) and for the guided one. The comparison of the flight paths is presented in Figure 5. In Figure 6, the errors between the target location and the actual position of the HSUAV are shown for the terminal phase of flight. Angular rates (roll, pitch, and yaw) are shown in Figure 7 and orientation angles in Figure 8.
Three-dimensional flight trajectory
Position errors (in the last phase of flight)
Angular rates
Euler angles (roll, pitch, and yaw)
The uncontrolled HSUAV ground range is about 8950.3 m, and the maximum achieved altitude is approximately 2680 m. For the unguided platform, the impact point coordinates are (6378.1, 6279.1) m, and the achieved miss distance was 77.58 m.
For the controlled one, the impact point is located at (6429.0, 6220.6) m, so only 2.18 m from the target. That means the lateral motors can successfully steer the HSUAV to the target.
All three components of the position error should ideally be 0. For the unguided object, the x and y errors (in the horizontal plane) deviated significantly from the desired value. On the other hand, for the guided one they are very small (<2 m).
The absolute value of the roll rate increases in the active portion of the flight (when the main motor operates). The minus sign in the first plot in Figure 7means that the HSUAV rotated counterclockwise when looking from aft. For the controlled flight, the oscillations of pitch and yaw rate signals after 30 s result from the firings of lateral motors.
For the guided object, the disturbances can also be observed on pitch and yaw angle time histories.
The forces generated by the lateral motors in the aeroballistic frame attached to the platform center of mass are presented in Figure 9. It can be observed that 21 lateral motors were used.
Forces generated by lateral motors
Later, a set of Monte-Carlo simulations were performed (300 runs for each scenario). The comparison of impact point dispersion for unguided and guided scenarios is presented in Figures 10 and 11.
Impact point dispersion (uncontrolled flights)
Impact point dispersion (controlled flights)
The CEP50 centered on the target for uncontrolled scenario one is 161.24 m. Some of the impact points are located 572.12 m from the aiming point (please see CEP100%).
There is a systematic offset between the mean point of impact (MPI) and the target. With the control system applied, CEP50% decreased to 4.50 m. That means using a control system reduces the impact point dispersion (measured using CEP50%) by approximately 35.83 times. The location of MPI coincides with the target. This is important from a practical point of view because the unintended collateral damage can be minimized.
Next, the C code was generated, and SIL simulations occurred. Finally, the obtained C code was implemented on the target hardware with ARM Cortex M4. The PIL simulations were evaluated, and the obtained results were compared with MIL and SIL. The comparison of outputs from the controller implemented in Simulink and using C code is presented in Figure 12. The abovementioned figure illustrates the number of firing commands obtained from the original model in Simulink and from the SIL and PIL simulations. In the ideal situation, lines should coincide with each other. Then, the original Simulink model and C code would be numerically equivalent. In practice, due to numerical rounding errors, the model and C code are numerically similar but not perfectly identical.
Number of activated lateral motors
Up to 30.37 s, no thruster was activated because the conditions given by Equation 28 were not fulfilled. That means the HSUAV was in an unguided flight. Then, a series of 21 firings occurred. The last motor in the sequence was activated in 41.01 s.
To better visualize and understand the eventual numerical errors, in Figure 13, the plot of time differences between Simulink and C code generated in SIL in evaluating firing commands is shown. It was assumed that the value of the acceptance threshold is 0.003 s.
Time differences in firing commands (Simulink vs C code)
The dots indicate that for most of the motors, there was no time difference between results from Simulink and C code. In 4 cases, the difference was smaller than 8e-15 s. The results obtained indicate that the C code possesses functionality identical to the original Simulink model.
In the next stage, a series of laboratory Hardware-in-the-loop experiments with the actual hardware were performed. This issue is not trivial because it is quite difficult to mimic the real conditions (for exam-ple, linear accelerations) in the static tests. Here, only the basic description of these tests is presented for legal reasons. The main goal of the tests was to ensure that the lateral motors could be fired in the right direction and at appropriate time intervals. For example, if the HSUAV should turn left (when looking from aft), then the thrusters must be fired on the right side of the fuselage. To verify the control unit, a simplified laboratory test stand was developed. The apparatus (Figure 14) consisted of a mounting, electric rotor, encoders, control unit, and operator’s control desk.
Control system prototype during the laboratory experiments
The main computer electronics were powered by the onboard rechargeable lithium polymer battery (12 V) integrated with the control system. The electric motor that drives the stand was connected to a different power source. The remaining elements (main solid propellant rocket motor and nose section) were removed because they are separate units (not neces-sary in the experiments). The developed control algorithm was implemented in the onboard computer.
The electric motor about the longitudinal axis rotated the control module with the same angular rate as obtained from the numerical simulation.
That way, the realistic input data about the roll angular rate and the roll angle from the navigation system were generated. The test stand is stationary, so the information about the position, linear speed, and accelerations cannot be delivered directly from the navigation module (values of these signals are zero). For that reason, the information about the other flight parameters was delivered directly to the onboard computer from the numerical simulation.
The firing commands were logged and compared with the data obtained from previous stages of testing (MIL, SIL). Firings of the lateral motors can be ineffective and problematic in laboratory conditions. A significant amount of time is required to prepare the test setup for the single test. After each test, the spent motors must be replaced with new ones. The lateral motors also produce a lot of smoke. To solve the abovementioned problems, instead of using real pulse thrusters, a set of LEDs (light-emitting diodes) was used to visualize the results (installed at the same locations as the original motors). The time of the operation of the single diode was set to be the same as the operation time of the real pulse thruster. The system was monitored using a high-speed camera (up to 5000 frames per second). Using this approach, the experiments can be repeated many times. That way, it was confirmed that the developed control system operates properly and is ready for flight trials. Real tests are required to prove the correctness of the system’s operation because even a sophisticated model cannot predict all the physical phenomena that can occur in flight.
Developing the software for the autopilot is complicated, challenging, and time-consuming. This process can be performed much more efficiently using a model-based design approach than traditional hand coding. Also, the integration and testing phases can be done in a structured and repeatable manner. The probability of introducing human-made errors can be minimized. It might be expected that this approach will gain more and more popularity in such applications. The main contribution of the reported research is the detailed description of model-based development software for the gasodynammically controlled guided HSUAV. Built-in MATLAB tools and customdeveloped verification procedures were used to create reliable software for the autopilot. SAVOIR guidelines were applied to perform the overall workflow. MIL, SIL, PIL, and HIL tests were completed to verify the correctness of the solution. The need for manual code implementation on the target hardware was minimized in the described process. Using MBD allowed moving some testing activities to the earlier stages of the system development. The results indicate that a good numerical similarity between the Simulink prototype and the C code was achieved. The presented results partially fill the existing literature gap and extend the simulation study reported in the papers [30,42].
Future research might involve code optimization to achieve higher computational efficiency and more experiments in laboratory conditions with real equipment. Hardware-in-the-Loop testing using dedicated Speedgoat target hardware and Simulink Real-Time might be evaluated. The performance of the developed control system will be carefully checked during the real flight tests.
The coordinate systems used in the mathematical model are presented in Figure 15.
Coordinate systems
The dynamic equations of motion were derived using the linear and angular momentum change the-orems for the 6 degrees of freedom rigid body taking into account the mass variation in time. In the body-fixed, non-inertial coordinate system Obxbybzb (the origin Ob does not coincide with the center of mass of the vehicle), the equations are as follows:
The nonlinear dynamic equations of motion of the projectile are as follows:
where 0 - zero matrix and 1 - unit matrix. In the short form, this is:
Quaternions were used to describe the object orientation:
The transformation matrix (39) could be used to formulate the relations between the rate of change of the position in Onxnynzn and linear velocities in Obxbybzb:
Quaternions could be converted to orientation angles (roll, pitch, and yaw) using the relations:
The initial quaternion could be calculated from the initial orientation angles in the following way:
The external forces and torques were calculated as the sum of aerodynamics
The vectors of aerodynamic force and moment are calculated with respect to point Oe, so with respect to point Ob they are given as:
Gravitational acceleration in the Onxnynzn coordinate system is given as
Thrust vector can deviate from the geometric lon-gitudinal axis of the HSUAV by angle ΘT in pitch plane and ΨT in the yaw plane, so the main motor loads are:
The gasodynamic control system is based on a set of identical correction lateral thrusters at the circum-ference of the body. The thrust and torque generated by the i-th thruster are:
The instantaneous mass of the vehicle is calculated as: