Linearized State-Space Model-Based Attitude Control for Rocket With Four Controllable Fins—Part 1-1: Basic Modeling and Identification

In this study, we present, analyze, and validate a linear state-space model for a simplified rocket with four fins. The model is validated experimentally using a wind tunnel. The proposed model consists of two main elements: a rocket body and fins. The body is represented as a cylinder integrated from infinitesimal sections along the longitudinal direction (x-axis), and each fin is interpreted as an individual rigid body. Applying a linear state space facilitates the use of linear analysis tools, leading to expedient system stability assessment and controller design. Stability analysis of the specific model is performed for validation. The results of the analysis showed that the phase margins of pitching and yawing were 16.9°, and that of rolling was 2.21°. However, the linearized state space carries the risk of discrepancies with the actual dynamics. The model is experimentally validated in a wind tunnel to reduce risk. The rocket model is affixed to a stationary jig in a wind tunnel and has four degrees of freedom. This arrangement allows for the verification of a coupled three-axis rotational model. The results show convergence within an average similarity of 77% in the linear range, confirming the reliability of the model. The validated model can be used for comprehensive analyses including control system design, performance optimization, and robustness analysis of a rocket with four fins in the future.


C d0,fi
Linear function offset of the i th fin's drag coefficient Total, aerodynamic, gravity force vector of the body N F bj Total force vector of the j th body element N F fi Total force vector of the i th fin N F g , F g,bj , F g,fi Gravitational force vector for the general system, j th element, i th fin, respectively N F l , F l,bj , F l,fi Lift force vector for the general system, j th element, i th fin, respectively Aerodynamical moment vector for the general system, j th element, i th fin, respectively With recent advancements in aerospace engineering, rockets designed for single use have transitioned to systems with diverse control objectives, including maneuver control, vertical takeoff, landing, and recovery [1], [2], [3].These systems require complex control inputs than traditional rocket controllers, which use the aerodynamic forces of the wings, attitude control thrusters, and thrust vectors of the main engines [4], [5], [6].These complex systems involve numerous states.A matrix-based state space is used to analyze and control these states [7], [8], [9].In this context, the use of a linearized state space is vital for understanding the dynamic characteristics of a system, facilitating stability, control system design, robustness, performance optimization, and transfer function analyses based on time-frequency analysis [10], [11].This study considers a linearized state-space model-based attitude control for a rocket that can be applied to larger systems with complex control inputs in future studies (models can be created, including thrust bias inputs or inputs such as canard fins and drag fins).
Linearization can encounter difficulties in models with complex subsystems [12].Over the years, multiple approaches for modeling the dynamics of rockets and missiles have been studied [13].
In this study, we first identified each component individually and detected the induced rotations during each experiment using a gimbal capable of a three-axis rotation in a wind tunnel.This allowed for the validation analysis of the proposed rocket system model.

A. MODELING
Papp [14] of the National University of Public Service at Budapest developed a mathematical model and system design for missiles and defined the flight conditions of air-to-air missiles.Papp designed mathematical models and controllers that matched specific flight conditions using fundamental aerodynamic coefficients such as lift and drag coefficients and trigonometric functions to create nonlinear systems for controller design.Although the mathematical design is relatively simple, it is yet to be validated.
Farhan [15] of the Weapons System Division of the Defense Science and Technology Organization of Australia reported a state-space model for autopilot design in the aerospace industry.A partially linearized nonlinear missile model was introduced.They used four cruciform wings arranged as control inputs and included subsystem models, including gyros and actuation servos.Although this study provides valuable insight for future controller designs, the complexity of the model makes it challenging to experimentally validate.
Biertumpfel et al. [16] of the Dresden University of Technology attempted a robustness analysis for launch vehicles ascending in the atmosphere.They presented a finite-horizon linear time-varying system assuming worst-case scenarios.The launch vehicle was controlled using thrust vectoring, and the rocket motion was described based on a nominal/reference trajectory.This study acknowledges the challenge of considering the effects of thrust variation.
Kisabo et al. [17], of the National Space Research and Development Agency at Suleja, decoupled the coupled nonlinear dynamics of a 6-DOF system and linearized it in state-space using the derivatives of the aerodynamic coefficients.The rocket model assumes thrust-vector control and uses four fins.The proposed model can be applied almost universally to aerospace vehicles by modifying the numerical values of the aerodynamic coefficients and their derivatives.
These models may exhibit discrepancies from actual flight dynamics because they have not been validated through actual flight tests and simulations.To bridge this gap, our study introduces the derivation of a rocket model with four controllable fins.We validated the model using a wind tunnel with hardware-in-the-loop simulation (HILS) [18] to establish a concrete testing methodology.

B. IDENTIFICATION
Miedzinski et al. [19] of the Warsaw University of Technology introduced a method for system identification using actual flight data from a suborbital rocket.Flight tests were conducted with biased inputs from the canards for test validation, and the test results were analyzed to validate the match between the model and actual data.One advantage of this method is that it is more accurate than verification using computational fluid dynamics.However, actual flight tests can be expensive and pose safety risks, necessitating system identification in a room.
Hann et al. [20] of the Canterbury University studied a sounding rocket in a vertical wind tunnel.In the experiment, rocket liftoff was simulated using a vertically erected wind tunnel, and the forward part was connected to a rotatable metal rod.This method effectively predicts the roll rate and angle by decoupling the disturbances from the intrinsic roll dynamics of the rocket frame.They used an integral-based parameter identification approach to consider wind disturbances equivalent to the movement of the actuator fins.This method is robust, computationally efficient, and accurately reflects the randomness of turbulent wind flow.However, the study only identified a single-axis system for roll control; the induced pitch and yaw rotations due to partial fin loss during the three-axis rotation in the wind tunnel were not confirmed.
Strub et al. [21] of the French-German Research Institute of Saint-Louis identified the pitch-axis rotation of a rocket-shaped projectile with controlled fins.Unlike previous roll axis rotation experiments, a rotating steel rod was connected to the center of gravity of the projectile in a horizontal wind tunnel, enabling pitch rotation.In this study, methods for installing a three-axis rotatable gimbal, collecting rotation angle data through an inertial measurement unit (IMU), controlling the bias angle of the fin with a servomotor, conducting experiments, and performing validation were evaluated.

C. CONTRIBUTION
Our study significantly contributes to rocket dynamics research by investigating, linearizing, deriving, and experimentally validating a model with fewer parameters representing rocket geometry.This study provides a solid foundation for future research and a practical approach to model validation and performance optimization through linearized state-space models.

D. PAPER ORGANIZATION
The remainder of this paper is organized as follows.Section II presents the definitions and derivation of the rocket model.Section III describes the model identification for simulation and experiments.Section IV presents the results of the stability analysis.Section V discusses the experiments and results.The conclusions are presented in Section VI.

II. MODELING
The bold symbols denote the matrices or vectorial quantities; italics denote the scalar quantities.
Rockets are typically designed to quickly exit high-density layers of the atmosphere at high acceleration rates to mitigate gravity and drag losses [22].However, owing to the inherent structural constraints of a rocket, the dynamic pressure applied to the rocket can only reach the designed maximum.Based on this premise, the proposed model assumes that the rocket maintains constant dynamic pressure throughout a specified flight trajectory.
For a Rocket with Four Controllable Fins, it is assumed that the rocket maintains a constant dynamic pressure after launch and ascends solely through fin control without propulsion.The rocket dynamics are illustrated in Fig. 1, where x, y, and z denote the linear displacements of the inertial reference frame (m), φ, θ, and ψ denote the rotation for the trajectory-relative frame (rad), and δ is the deflection of the fin as input (rad).Each fin is installed at position L f , which is the location vector from the body frame.The origin of the body frame was the center of mass (CM).
The launch site was considered the inertial reference frame.The rocket progressed along the trajectory based on the rela-tive wind; thus, x −x 0 represents the distance from the launch site along the trajectory.

A. SYSTEM MODEL 1) RIGID BODY MODEL
The trajectory-relative motion vector r and trajectory-relative attitude vectors can be expressed as the sum of each component: where the m represents the total mass (kg) of the rocket, F a,fi and F g,fi represent the i th fin's aerodynamic and gravity force vectors (N ), respectively, F a,bj and F a,bj represent the j th element's aerodynamic and gravity force vectors (N ), respectively, M a,fi and M g,fi represent the i th fin's aerodynamic and gravity moment vectors (N • m), M a,bj and M g,bj represent the j th element's aerodynamic and gravity moment vectors (N • m), and J represents the total inertia tensor (kg about the CM of the rocket, expressed as In particular, in Eq. ( 1) and ( 2), i represents the set of installed fins, and j represents the set of elements for the Riemann integral representation of the body of the rocket 146018 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
from the front L x,m to the end location L x,m − L h,b , expressed as the M g,fi and M g,bj can be neglected when the rocket is assumed to be a rigid body, and gravity is a constant CM, defined as the center of gravity, expressed as where M g is the total moment of the rocket, defined as 2) AERODYNAMIC MODEL The drag coefficient C d , sideslip coefficient C s , and lift coefficient C l are defined as follows [23]: where F d , F s , and F l represent the axial drag, sideslip, and lift forces (N ), respectively, A x , A y , and A z represent the reference cross-sectional areas (m 2 ) perpendicular to the x b , y b , and z b axes, respectively.The dynamic pressure q is defined as [24]: where ρ is the air density (kg • m −3 ), and v is the air speed (m • s −1 ).The aerodynamic coefficients related to the effective angle of attack α and effective sideslip angle β can typically be linearized to partial derivatives [25], [26].Thus, using Eq. ( 8), the general equations for the drag, lift, and sideslip forces can be represented by the effective angle of attack α and sideslip angle β as given in Eq. (10).
where C dα and C dβ represent the linearized aerodynamic partial derivatives of the drag coefficient for α and β (No Unit), respectively, C lα and C sβ denote the partial derivatives of the lift coefficient for α and β (No Unit), respectively, C d0 , C l0 , and C s0 represent the linear function offsets of the drag, lift, and sideslip force coefficients (No Unit), respectively.Aerodynamic forces F a,bj and F a,fi for j and i, respectively, are defined as where the rotation matrix denoted as R b f fi converts the i th fin's local frame to the body frame using a direction cosine matrix (DCM).Rotation matrices for fin are designed as follows: where L fi denotes the location of the center of pressure (CP) in the fin.Thus, R b f fi is calculated using the DCM as The aerodynamic moment can be expressed as the cross-product of the force and displacement between the CM and the point of action of the force as follows: where the location vectors of the CP for j and i are denoted as L bj and L fi and each CP acts as a point action of forces F a,bj and F a,fi , respectively.

B. STATE-SPACE MATRIX DERIVATION
The total force of the body F b and moment M b were derived using Eq. ( 15) and ( 16) as follows: where F b and M b are the total force and moment of the body, respectively, and F g,bj and M g,bj are the gravitational force and moment, respectively.Since a gimbal is attached to the CM of each rocket, Thus, the total force and moment of the fins F f and M f are summarized as follows: The mass of the fins is considered to be included in the body; thus, the gravitational force and moment caused by the fins are neglected.
Because the rocket is axisymmetric for the x b axis, the coefficients of the lift and sideslip linear functions C dα,bj and C lα,bj are identical, and the linear function offset can be neglected.
The rocket fins were assumed to use a symmetric airfoil.Therefore, the sideslip coefficient and linear function offset of the thin airfoil were neglected.
The coefficients C dα,bj , C dβ,bj , C dα,fi , and C dβ,fi of the linear function representing the drag coefficient were ignored according to the cosine approximation method.
As shown in Fig. 1, Eq. ( 25), each i th fin is installed axisymmetrically about x b , and L x,fi is assumed to be identical.
The body is an integration of infinitesimal sections, each with a thickness of L h,bj and location vector L bj , defined as The reference sectional areas (m 2 ) A x , A y , and A z for bodies A x,b , A y,b , and A z,b , are defined as Thus, L 1 and L 2 for the area multiplied by each element's location, L x,bj and square, respectively, are expressed as follows: The surface areas of fins S fi are all considered identical to S. Thus, A x , A y , and A z for fins A y,fi , A y,fi , and A z,fi , respectively, are defined as follows: The q is approximated by ẋ and can be derived as When the wind velocity of the inertial reference frame assumed zero, and the velocity ẋ is significantly greater than ż and ẏ, the effective angles of attack (α bj and α fi ) and effective sideslip angles (β bj and β fi ) for the j th body element and i th fin, respectively, can be determined using an inverse tangent approximation [27], as shown in Eq. (32). where, the v c is the constant wind speed for the linearization of ẋ; the velocity of the y w axes (ẏ bj and ẏfi ) and z w axes (ż bj and żfi ) for the j th body elements and i th fin, respectively, are summarized in Eq. ( 34).The deflect angle of the i th fin is δ i and changes the designed angle of attack.δ i is considered as an element set of the input U.Each α f i is summarized in Eq. ( 35): Thus F b , F f , M b , and M f can be summarized as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

M
The detailed derivation process is presented in Appendix A.
In the state-space representation, r and can be summarized as Using Eq. ( 36), (37), ( 38), ( 39), ( 40), (41), the state-space matrix is derived as shown in Eq. ( 42), as shown at the bottom of the page, where, as (43), shown at the bottom of the page, where A 22 and A 42 are related to ṙ, and according to Eq. ( 20) were considered negligible. (44)

III. IDENTIFICATION
Wind tunnel experiments were conducted to obtain the coefficients, and the linearized partial derivatives of the lift and drag were calculated for each angle of attack, as shown in Fig. 2. Fig. 6(a) in Section V shows the identification for C lα,fi , C lα,bj , C sβ,bj , and C d0,bj are identified in the configuration as shown in Fig. 6(b).Each coefficient in the result is a linear function that includes the derivative coefficients.The linearization scope is selected from 0 • to 30 • or 0.52 radians for the linear approximation of trigonometric functions.We selected the linearized lift coefficient C lα,fi = 3.261 (rad −1 ) from the wind tunnel results; C lα,bj = C sβ,bj = 0.726 (rad −1 ), C d0,bj = 0.5.The values of S fi , L h,b , L w,b , L fi , and L m were measured directly.We refer to the nomenclature for all values.

IV. STABILITY ANALYSIS A. TIME RESPONSE
Stability analysis was conducted for each axis using a linear state space.The rocket was axisymmetric, yielding similar results for the pitch and yaw in both axes.Each step input used U θ , U ψ , and U φ , with corresponding outputs C θ , C ψ , and C φ , defined in Eq. ( 45) for the pitch, yaw, and roll rotations, respectively.This approach prevents unintended coupled rotations owing to asymmetric wing deflections.The input and output U and C in Eq. ( 42) were replaced with the following in the analysis: The pitch and yaw step responses are shown in Fig. 3(a).Fig. 3(b) shows significantly different results for the roll.For the pitch and yaw responses, there is convergence of the output for unit input 1; however, there is divergence in the roll.Furthermore, the outputs for the pitch and yaw converged in the negative direction of the input.Based on these findings, the input inversion may be necessary when analyzing Bode plots.

B. BODE PLOT
A Bode plot graphically represents the frequency-response characteristics and can be used to analyze and test the stability of the feedback control system by obtaining the gain margin (GM), phase margin (PM), and delay margin (DM) of the amplifier.Although the plant in this study does not include a controller, the stability analysis of the system allows for the evaluation of the classic wing stabilization control of the rocket and the determination of the required controller performance.
The PM is the difference between the phase shift and -180 • ; the gain is 0 dB in the Bode plot.GM is the difference between the gain and 0 dB when the phase shift is -180 • .
As shown in Fig. 4(a), the PM of the system for the pitch and yaw rotations was 16.9 • at 3.87 rad/s.GM was 0.807 dB at a phase crossover frequency of 180 • .For the roll rotation, as illustrated in Fig. 4(b), PM was 2.21 • at 4.63 rad/s, and no GM was observed.
It is generally known that if the PM is less than 30 • , the system becomes unstable and reacts significantly.If the angle was greater than 60 • , the system was likely to have a slower response speed.In general, a PM angle of 45 • is considered appropriate [28].In this system, for the pitch and yaw, the PM was 16.9 • when the magnitude was 0 dB.Although the results converged, as shown in Fig. 4(a), the system was unstable and required additional stability or controller.For roll rotation, Fig. 4(b) shows the unstable system with a PM of 2.21 • , significantly diverging in the 30 • -60 • range owing to the input.The DM for pitch and yaw were analyzed as 76.8 ms and 8.34 ms for roll.

C. NYQUIST PLOT
The Nyquist criterion can be used in cases where the transfer function of a control system includes error factors arising from experimental methods, to confirm the relative stability or intuitively understand the influence of the characteristics of a particular control element on the stability of the entire system.
According to the Nyquist criterion, Fig. 5(a) shows the stability as the mapping rotates clockwise and is positioned to the right of pole (-1,0) [29].In contrast, Fig. 5(b) shows instability when the mapping is located to the left of the pole (0,0).Based on this criterion, although the predicted system for the pitch and yaw inherently exhibited stability, the roll was unstable.

V. EXPERIMENT AND RESULTS
The experimental setup is shown in Fig. 6.The rocket model was mounted on a 3-DOF rotatable gimbal, equipped with a servo motor for adjusting δ i , a computer, and an IMU for collecting attitude information.Parameter identification was performed using the load cells shown in Fig. 6(a) and Fig. 6 (b) and a series of experimental processes, and the results were applied to the model.

A. EXPERIMENT
The gimbal was designed to minimize mass and inertia effects without additional power and was installed at the center of gravity for rocket rotation.As a design limitation, a bearing larger than the rocket was used for roll rotation, possibly leading to a higher friction than the smaller bearings used for pitch and yaw rotation.A nacelle was installed to prevent aerodynamic effects of the gimbal on the wings.
Computational fluid dynamics and actual measurements confirmed consistent wind-speed delivery to the fin.However, subtle vibrations can occur owing to disturbances and rotational angle limitations, which can affect the results.
To minimize these effects, the peak angles were analyzed for the pitch and yaw to ensure that they exceeded the angle limitations caused by the nacelle, as presented in Section IV (stability analysis), and the inputs were adjusted accordingly.For roll rotation, the step response continued to diverge; φ reached 1200 • /s at a 10 • average δ i .
However, it was experimentally confirmed that the rotation stopped if the deflection angle decreased to less than 15 • , owing to limitations in the experimental setup, such as the gimbal static friction.Conversely, the installed IMU could not measure rotations exceeding 2000 • /s; thus, appropriately adjusted values were used in the validation test.The input The input U (t) was fed simultaneously into the model for software-in-the-loop simulation (SILS) and HILS devices in real-time, as shown in Fig. 6.The errors e (t) were measured by comparing the outputs Ŷ (t) and Y (t) from the SILS and HILS, respectively.
where Ŷ and Y were calculated using Eq.(42).However, the output matrix is replaced by Eq. (47).Consequently, the output is expressed as a 1 × 1 scalar matrix, and the error e according to the roll, pitch, and yaw is scalar.

B. RESULTS
Each experiment for the input scenarios in Eq. ( 47) is shown in Fig. 7.The SILS and HILS plots are compared in Fig. 7(a), (c), and (e).The ERR was calculated using Eq.(48) in  as shown in the analysis results in Fig. 3(c) and the experimental results in 7(e); thus, it was not possible to calculate the maximum and steady-state responses according to the input.Instead, the maximum values and amplitudes of the angular velocity were calculated, as shown in Fig. 3(d) and 7(g).The CMD of each experiment was not changed until 13 sec for the result to settle based on the settling time expected from the step response analysis results.Table 1 summarizes the peaks Ŷp and Y p , and the steady-state Ŷs and Y s of SILS and HILS, respectively, e p and e m denote the absolute maximum peak error and root mean square (RMS), respectively.Fit percentage low and For the roll, the fit percentage was calculated based on the angular velocity and errors.
In the experiment result, the Y p values of pitch and yaw were 12.5 • and 9.8 • , respectively, and the e p values were 4.9 • and 6.8 • , respectively.This result indicates that unlike the axisymmetric model, the Ŷp result is similar to the pitch and yaw, and the experimental rocket does not have the same perfect axisymmetry.The RMS errors were relatively similar to the peak errors at 1.4 • and 1.8 • for the pitch, yaw, and roll, and the calculated average fits were 86.3% and 76.9%, respectively.In contrast, the minimum fits were 61.2% and 30.4% for the peak standard pitch and yaw, respectively.For roll, it took approximately 0.6 s for the rocket to rotate 180 • ; the peak error based on 0-0.6 s was 37.7 • .The peak angular velocity was 1483 • /s from 0 to 30 s.The resulting error was 1343.6 • /s; the average fit was 68.0%, whereas the minimum fit was 9.4%.

VI. CONCLUSION AND FUTURE RESEARCH
Overall, the model similarity was 77%, however, the similarity for the peak error in some experiments was as low as 10%.Excluding the roll experiment, the average fit was 81.6% and the lowest fit was 30.4%.Excluding the yaw experiment, the lowest fit was approximately 61.2%.This was due to the limitations of the experimental configuration, which assumed that L w,b is constant during axisymmetry.In contrast, the body of the experimental device differed by 0.05 in the largest case and by 0.01 in the smallest case.However, it is difficult to solve a problem based only on this effect, and the aerodynamic effects of the gimbal and wind tunnel nacelles must be considered.In addition, for the roll, the frictional force on the rotating bearing can affect static and high speeds.
Nevertheless, this study confirmed the similarity of the model and the limitations of the model and experiments.Therefore, it would be beneficial to improve the model and conduct additional experiments.
This series is divided into four parts.Each part of the system is configured as shown in Appendix B.
Thus, F b , F f , M b , and M f are derived as where the cross-sections A x,bj , A y,bj , A z,bj , A x,fi , A y,fi , and A z,fi for each axis are defined as follows: Each cross-section is given by Eq. ( 27) and (30) as follows: Thus, Eq. ( 52), (53), (54), and (55) can be organized as 146026 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Because C lα,bj is a constant, it can be summarized as ( 71) and ( 72), shown at the bottom of the page.
In Eq. ( 28) and ( 29), L 1 and L 2 are defined as  78)and (79) can be summarized as follows: Eq. ( 35), from Eq. ( 33) is expressed as L x,fi is a constant, defined in Eq. (25) as Forces F f and moment M f of the fins are derived as

J 2 j
Total moment of inertia tensor kg•m 2 J xx Moment of inertia for rotation around the x-axis 0.023 kg•m 2 J xy Moment of inertia around the x-axis when the object rotates around the y-axis 0 kg•m 2 J xz Moment of inertia around the x-axis when the object rotates around the z-axis 0 kg•m 2 J yx Moment of inertia around the y-axis when the object rotates around the x-axis 0 kg•m 2 J yy Moment of inertia for rotation around the y-axis 0.436 kg•m 2 J yz Moment of inertia around the y-axis when the object rotates around the z-axis 0 kg•m 2 J zx Moment of inertia around the z-axis when the object rotates around the x-axis 0 kg•m 2 J zy Moment of inertia around the z-axis when the object rotates around the y-axis 0 kg•m 2 J zz Moment of inertia for rotation around the z-axis 0.436 kg•m Element number of the body (for j ∈ L m x , L m x − L b h ) No Unit L bj Position vector of the j th element relative to the body frame m L x,bj Position of the j th element on the x-axis m L y,bj Position of the j th element on the y-axis 0 mm L z,bj Position of the j th element on the z-axis 0 m L h,b Length of the rocket body 1.550 m L h,bj Length of the j th element of the body m L w,b Reference diameter of the rocket body 0.050 m L fi Position vector of the i th fin relative to the body frame m L x,fi Position of the i th fin on the xPosition of the body front on the x b -axis 0.860 m L y,m , L z,m Position of the body front on the y b -and z b -axis, respectively 0 m M a , M a,bj , M a,fi

FIGURE 1 .
FIGURE 1. Free-body diagram of simplified rocket with four fins.

FIGURE 2 .
FIGURE 2. Wind tunnel result of fin's lift coefficient by angle of attack.

FIGURE 6 .
FIGURE 6. Validation process via experiment: (a) identification for fin; (b) identification for body and experiment for validation.

Fig. 7 (
b), (c), (d).Each figure shows the input command CMD as an applied maximum δ fi value in Eq. (47).Fig. 7(e) and (f) show graphs comparing the angular accelerations and their errors in the roll rotation.The roll rotation diverged,