Hardware In-the-Loop Simulation for Control System

Aerospace Indonesia Meeting, Bandung, Indonesia 27 July 2005

1

Hardware In-the-Loop Simulation for Control System Designs of Model Helicopter T. Sudiyanto† ,A. Budiyono‡ and H.Y. Sutarto§ Abstract– Many control system designs have been developed for controlling the dynamic process of aerospace vehicles. Yet, among those designs, few are readily suitable for implementation. Moreover, most of the applicable aerospace control designs are built through flight tests that perform real processes on available real vehicles. The need to use real vehicles in the control design cycle poses a high risk and is prohibitively costly. To solve this problem, a real-time simulation concept that employs cheap, practical, and rapid-to-build modular hardware, which can simulate a nearly real process in lab’s environment, was proposed in this research work. The paper elaborated the integration of a nonlinear model of a small-scale helicopter, Yamaha R-50, with the sensor and actuator models into a modular Simulink model. The model was then uploaded into a dSPACE’s digital signal processor (DSP) chip. Based on that Simulink model, the DSP performed a real-time computing process and behaved like the represented helicopter dynamic system. The overall scheme of Hardware-In-the-Loop simulation system can play a significant role in the efficient control system design for aerial vehicles. Keywords– Hardware-In-the-Loop simulation, sensor model, Kalman filter, helicopter dynamics modeling

A1 , B1 a1 , b1

a1s , b1s ab B bb cb CD0 d e Fax , Fnorm Fx , Fy , Fz fwake G(s) h Ib

Ix , Iy , Iz , Ixz is k1 , k2 , kc k, l kCR/SP kMR/SP ,kMR/β

Mx , My , Mz mg P p, q, r Q R Sxeff , Syeff , Szeff T u, v, w (u, v, w)B (u, v, w)H

Nomenclature lateral and longitudinal rotor blade feathering coefficients longitudinal and lateral rotor blade flapping coefficients relative to swash-plate plane longitudinal and lateral rotor blade flapping coefficients relative to hub plane blade lift slope number of blades blades’ plane span length blade’s chord length effective blade drag coefficients horizontal distance to the center of gravity rotor hinge offset axial and normal aerodynamic force force components in body frame wake function for low/high speed transfer function in Laplace domain vertical distance to the center of gravity rotor blade’s moment of inertia

† Department of Aeronautics and Astronautics, Bandung Institute of Technology, Jl. Ganesha 10, Bandung 40132, Indonesia [email protected] ‡ Department of Aeronautics and Astronautics, Bandung Institute of Technology, Jl. Ganesha 10, Bandung 40132, Indonesia [email protected]. The author to whom all correspondce should be addressed. § Department of Aeronautics and Astronautics, Bandung Institute of Technology, Jl. Ganesha 10, Bandung 40132, Indonesia [email protected]

V (v, w)i (v, w)b (v, w)r α βC, βS γ δ3 θ0 λ ξ ρ φ, θ, ψ Ω Ωf ω in , ω off

rigid body’s moment of inertia rotor shaft tilt cross coupling coefficients directional parameters for lift generated by tail planes coefficient defining control rotor blade pitch due to swash-plate tilt coefficient defining main rotor blade pitch due to swash-plate tilt and due to control rotor tilt moment components in body frame helicopter’s weight rotor power angular rates in body frame rotor torque rotor radius effective wet area in body frame rotor thrust velocity components relative to air in body frame velocity components in body frame velocity components in local horizon frame total air speed induced velocity due to tail rotor and main rotor thrust tail rotor and main rotor average velocity relative to air tail rotor and main rotor disk velocity relative to air tail plane’s local angle of attack longitudinal and lateral control rotor tilt Lock number delta-three-angle rotor blade pitch constant directional parameter for rotating rotor limited extension parameter air density attitude angle rotor rotational speed coefficient defining change in main rotor natural frequency due to hinge offset flap rate coefficients for in-axis and offaxis motion I. Introduction

A

N Implementation and testing of control systems by a hardware-in-the-loop (HIL) simulation in aerospace

2

engineering is increasingly being required for the design as it becomes a very beneficial tool in acquiring ’real’ data without taking a risk of loosing any real vehicles. Another advantage of this simulation is that theoretically for a successful controller test, the mathematical helicopter model simply need to be replaced by the real helicopter. HIL simulation is characterized by the operation of real components in connection with real-time simulated components. Usually, the control system hardware and software is the real system while the controlled plant can be either fully or partially simulated. The controlled plant discussed in this paper consists of actuators, a nonlinear model of a small-scale helicopter, Yamaha R-50, and sensors. The control system is excluded from this paper’s context. The accuracy of the mathematical model of the helicopter dynamic becomes the major concern in the HIL simulation. The better the mathematical model of the whole system and all important aspects such as disturbances, noise, and time delays, the easier is the transition from simulation to real flight, and more important, the lower the risk of loosing the aircraft because of system failure. The mathematical model of the Yamaha R-50 helicopter used in this work is already derived and available in [1]. The GPS model is available in [2].

       Mx Mx Mx Mx  My  =  My  +  My  +  My  + Mz Mz MR Mz TR Mz fus     Mx Mx  My  +  My  Mz HT Mz VT These equations are in general valid over the entire flight envelope for the typical aircraft. 

A. Main Rotor Model The following are the force and moment equations for the main rotor: FxMR = −TMR sin (a1s + is )

(10)

FyMR = TMR sin b1s

(11)

FzMR = −TMR cos (a1s + is ) cos b1s

(12)

MxMR

= FyMR hMR +

II. Helicopter Mathematical Model

λ

For manned and unmanned simulation of flight, the six rigid body DOF is accurate in describing the system response for flight mechanical and controlling purposes. Therefore, the helicopter model is based on the equations of unsteady motion of the rigid body derived in body axis system. ³ · ´ Fx = (mg sin θ) + m uB + qwB − rvB ³ · ´ X Fy = − (mg sin φ cos θ) + m vB − pwB + ruB ³ · ´ X Fz = − (mg cos φ cos θ) + m wB + pvB − quB ´ ³· X · Mx = Ix p + (Iz − Iy ) qr − Ixz r + pq X ¡ ¢ · My = Iy q + (Ix − Iz ) pr − Ixz p2 − r2 ³· ´ X · Mz = Iz r + (Iy − Ix ) pq − Ixz p − qr X

·

φ = p + (q sin φ + r cos φ) tan θ; θ 6= ± ·

π 2

λ

dMyMR a1s + da1

dM (−b1s + B1MR − k1 a1s ) dB1

(14)

(2)

eMR where k1 = tan δ 3 ≡ pitch horn length According to the above equations the three force components of the main rotor comes from the thrust generated by the rotating blades. The roll and pitch moments are generated by the thrust, flapping angle with respect to main rotor hub plane and feathering angle of the rotating blades. The yaw moment is contributed only by main rotor torque. The blades flapping angle coefficients with respect to hub plane equations in term of the blade flapping angle coefficients with respect to the swash-plate plane are given in the form

(3) (4) (5) (6) (7)

ψ = (q sin φ + r cos φ) sec θ

(9)



= −Fx hMR + FzMR dMR +

(13)

MzMR = λQMR

(8)

where   Fx  Fy  = Fz    Fx  Fy  +  Fz HT and

dMxMR (a1s + A1MR − k1 b1s ) dA1

(1)

θ = q cos φ − r sin φ ·

MyMR

dMxMR b1s + db1

     Fx Fx Fx  Fy  +  Fy  +  Fy  + Fz MR Fz TR Fz fus  Fx Fy  Fz VT

a1s = a1 − B1MR + λkc b1s +

b1s = b1 + A1MR − λkc a1s +

da1s u (1 + 2fwake ) du

(16)

db1s v (1 + fwake ) dv

(17)

ΩMR where kc = k1 + k2 , k2 = 34 ReMR Ωf , MR ³ ´ γ ΩMR and Ωf = MR16 1 + 83 ReMR , MR ³ ´ ρab cb R where γ MR = Ib MR

(15)

3

The equations of main rotor flapping motion relative to the swash-plate plane is in the form ·

a1 = −ω in a1 − λω off b1 − q ·

b1 = −ω in b1 − λω off a1 − p where ω off =

³ΩMR ´ 2 ΩMR 1+ Ω f ΩMR ω off Ωf

(18)

wr = (a1s + is ) u − b1s v + w V 2 = u2 + v2 + w2 (wr − 2wi ) The main rotor torque is computed by

(19)

,

QMR =

·

where a1 = 0 and b1 = 0 can be used. Solving (16) and (17) for a1 and b1 , equation (18) (19) can be solved. The longitudinal and lateral swash-plate tilt B1SP and A1SP commanded by the pilot is distributed by mechanical linkage to the control rotor as well as directly to the main rotor. The resulting tilt of the control rotor plane due to control blade pitch changes commands the main rotor blade pitch through mechanical linkage. Then, the resulting main rotor input is defined as

PMR = Ppr + Pi + Pc + Ppa

B1MR = kMR/SP B1SP − kMR/β β CCR

(20)

A1MR = kMR/SP A1SP + kMR/β β SCR

(21)

The computed moment derivatives due to blade feathering and blade flapping is given as

=

dMxMR db1s

¶ µ 1 eMR dMyMR 1 = B × dB1 2 6 RMR ³ ´ 2 2 ρab cb RMR (ΩMR RMR )

¶ µ dMyMR 3 eMR = = Bλ × da1s 4 RMR Ã ! 2 2 ρab cb RMR (ΩMR RMR ) γ MR

(22)

(27)

where Ppr

=

µ

¶ 1 1 ρ CD0 bcb RMR (ΩMR RMR ) × 2 4 h ¡ ¢i 2 (ΩMR RMR ) + 4.6 u2 + v2 Pi = TMR wi ·

dMxMR dA1

(26)

The main rotor power

and ω in = Compared to the rigid body dynamic, the main rotor dynamic are usually very fast. Therefore the conditions ·

PMR ΩMR

Pc = mg z H ³ ´ Ppa = − Fxfus u + Fyfus v + Fzfus (w − wi )

(28) (29) (30)

(31)

Control Rotor Model Since the control rotor (some references refer this component as stabilizer bar) has only small profile at the end of the bar, its contribution to the helicopter resultant force is negligible, while the aerodynamic moment cause the control rotor plane to tilt. The final equations representing the control rotor dynamic are ·

1 ΩCR γ CR ξ × 16 ¶ q − A1CR + =0 ΩCR

β SCR − p +

µ β SCR + ∆β SCR

(23)

(32)

where The main rotor thrust is computed by solving these equations iteratively on every time step: 1 2 TMR = (wb − wi ) · B · ρΩMR RMR ab cb 4

wi2

=

sµ

1 2 V 2

¶2

+

µ

1 TMR 2 ) 2 ρ (πRMR

¶2

1 − V2 2

∆β CCR =

θ0CR u γ CR ξΩCR RCR

(33)

∆β SCR =

θ0CR v γ CR ξΩCR RCR

(34)

(24)

(25)

are the additional flapping angle due to relative wind velocity. where

where µ ¶ 2 3 wb = wr + ΩMR RMR θColl + θtwist 3 4

µ ¶4 lb ξ =1− 1− RCR

(35)

4

And the equations of control rotor pitch angle as used in (20) and (21) are B1CR = kCR/SP B1CR

B. Fuselage Model The force and moment equations for the fuselage is computed as

(36)

A1CR = kCR/SP A1CR

³ ´ 1 Fxfus = − ρu2 Sxeff 2 fus

(49)

³ ´ 1 Fyfus = − ρv 2 Syeff 2 fus

(50)

³ ´ 1 Fzfus = − ρw2 Szeff 2 fus

(51)

Mxfus = Fyfus hcpfus

(52)

Myfus = −Fxfus hcpfus + Fzfus dcpfus

(53)

Mzfus = −Fyfus dcpfus

(54)

(37)

Tail Rotor Model The force and moment equations for the tail rotor is acquired by the same manner as they are for the main rotor, except that the tail rotor blades are not hinged so they can not flap and that it can only be commanded by deflecting its blades collective pitch angle. FxTR FyTR FzTR MxTR MyTR MzTR

= 0 = −TTR = 0

(38) (39) (40)

= FyTR hTR = λQTR = −FyTR dTR

(41) (42) (43)

C. Horizontal Tail Plane

The thrust is derived iteratively from these equations

vi2

=

1 2 TTR = (vb − vi ) · B · ρΩTR RTR abTR cbTR 4

(44)

sµ

(45)

1 2 V 2 TR

¶2

+

µ

1 TTR 2 ) 2 ρ (σ TR πRTR

¶2

1 2 − VTR 2

The horizontal tail plane has a symmetric cross section, and only the normal force is considered since the axial force is negligible. The force and moment equations for the horizontal tail plane is computed as

FxHT FyHT Fxfus

where µ ¶ 2 3 vb = vr + ΩTR RTR θCollT R + θtwistTR 3 4

VT2R = u2 + v2 (vr − 2vi ) + w2

PprTR

(46)

µ

¶ 1 1 = ρ CD0TR bTR cbTR RTR (ΩTR RTR ) × 2 4 h ¡ ¢i 2 (ΩTR RTR ) + 4.6 u2 + w2 (47) Pi = TTR vi

= 0 = Fxfus dHT = 0

(58) (59) (60)

where

The tail rotor torque is computed by PTR ΩTR

(55) (56) (57)

MxHT MyHT MzHT

vr = v

QTR =

= FaxHT = 0 = 0 = FnormHT = Lif tHT cos αHT

(48)

¢ 1 ¡ 2 SeffHT ; Lif tHT = k · ρ u2HT + wHT 2 0 ≤ |αHT | < |αHTstall | ¢ 1 ¡ 2 Lif tHT = k · ρ u2HT + wHT (SeffHT )stall ; 2 |αHT | ≥ |αHTstall | k k

= −1; w > 0 = 1; w < 0

(61)

(62)

5

D. Vertical Tail Plane Like the horizontal tail, the vertical tail plane has a symmetric cross section, and only the normal force is considered since the axial force is negligible. The force and moment equations for the horizontal tail plane is computed as

FxVT FyVT FxVT

= FaxVT = 0 = FnormVT = Lif tVT cos αVT = 0

(63) (64) (65)

MxVT My

= FyVT hVT = 0

(66) (67)

Mz

= −Fyfus dVT

(68)

VT VT

where ¢ 1 ¡ 2 Lif tVT = l · ρ u2VT + wVT SeffVT ; 2 0 ≤ |αVT | < |αVTstall | ¢ 1 ¡ 2 SeffVT ; Lif tVT = l · ρ u2VT + wVT 2 0 ≤ |αVT | < |αVTstall |

(69)

(70)

E. GPS Model The Global Positioning System is a satellite navigation system that allows the user to acquire accurate determination of position and velocity based on noisy observation of the satellite signals. For the purpose of this simulation, a simple model has been developed to generate the errors of the GPS system and a filter to reduce the noise intensity of the GPS’ output. The models of both components are included in the position and velocity channels. The model has the same structure for both position and velocity but with different parameter values. The block diagram is presented in figure 5 and the values of the parameter of the model is in table 71. The main characteristic of the GPS which have been considered are latency, update rate, accuracy and error dynamics parameter. The update rate represents the rate at which the position and velocity signals are sent to the receiving processor and is modeled as quantization. The latency is the time delay that occurs between the time the satellite information is received and the time the position or velocity output is sent to the receiver. It is modelled as a pure time delay. The accuracy is the radius of the circle with the origin at the actual position or velocity which contains 50% of the sensors output values. The error of the GPS sensor package are generated as output of a first order linear differential equation with random Gaussian input and initial condition.

Update Rate, T2 Latency, τ Accuracy, k Error Dynamics Parameter, a

Position 5 Hz 0.075 s 0.65 ft 0.5 s

Velocity 5 Hz 0.075 s 0.1 ft 2.5 s (71)

A discrete Kalman Filter Algorithm is implemented to the GPS model. The algorithm is illustrated in figure 6 Detailed discussion about Kalman filtering can be observed in [3]. The definition of the algorithm’s parameters are as this followings. b− Pra-estimated state, x j bj Estimated state, x b− Next cycle pra-estimated state, x j+1 Error covariance of the pra-estimated state, P− j Error covariance of the estimated state, Pj Error covariance of the next cycle pra-estimated state, P− j+1 Kalman Gain, Kj Measured state, zj State transition coefficient, Φ Measurement coefficient, H Error covariance of the proccess noise, Q Error covariance of the measurement noise, R

F. Actuator Model Four identical actuators are used to control the four input channels. Two actuators control the swash-plate, and the other two control the main and tail rotor collective pitch. The actuators’ dynamic is approximated with a firstorder transfer function:

Gact(s) =

15 s + 15

(72)

III. The HIL Simulation A. Simulation Configuration The real-time HIL Simulation test facility is simply illustrated in figure 1. The upper block represents a plant on which the controlling device being tested. The mathematical model described in the previous section is then uploaded into this block and resides as its software as seen in figure 2. Figure 3 shows the elements of the model software. Each elements represents a helicopter’s component module. B. Simulation Procedure First of all, the helicopter mathematical model described in the previous section is built into a Simulink model. Once the model becomes executable in off-line simulation, it can be uploaded into a dSPACE chip. A complete sequence of this procedure is shown in figure 4. On the off-line

6

simulation, the values of the force and the moment resultants in z-axis is constrained to zero. Input values corresponding to this constraint are plotted in figure 7. The constant values as they reach a steady state will serve as the initial input set to run the off-line test of this model. The helicopter’s response these constant trim val   toward  θCollMR 6.11 ◦  θCollTR   11.38 ◦      ues   B1SP  =  0 ◦  are plotted in fig 8, 9, A1SP 0◦ 10, 11, 12, and 13. This off-line simulation results will be used as a comparator for preliminary analysis. IV. Discussion on Results The θCollMR has a direct effect in controlling the helicopter’s vertical motion. Hence, the z-axis component of its translational time response has the shortest time-to-steady compared to the other two components. Meanwhile, the main rotor rotation produces a yawing moment which need to be balanced. The balancing moment is generated by the θCollTR as it produces thrust at a horizontal distance from the helicopter center of gravity. Also, the tail rotor vertical distance makes it generate a rolling moment, which need to be balanced. Meanwhile, the tail rotor rotation produces a pitching moment which need to be balanced. This rolling and pitching moments provoke the main rotor tip path plane about its axis which counteracts those moments. The real-time simulation results with the dSPACE hardware set is not available yet due to the compatibility problem between the helicopter model built in Simulink 5 with the Simulink 2.2 on the dSPACE’s host PC. The complete report of this research including its real-time simulation results will soon be available.

Fig. 1. HIL Simulation Loop

V. Concluding Remarks All outputs except the heading angle were able to reach their steady state without any active control. At this point we may conclude that the modeled helicopter has a stable characteristic on at least one flight configuration. [5] References [1] C. Munzinger, “Development of a real-time flight simulator for an experimental model helicopter,” Master’s thesis, School of Aerospace Engineering, Georgia Institute of Technology, 1998. [2] M. G. Perhinschi and J. V. R. Prasad, “A simulation model of an autonomous helicopter,” [3] R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman Filtering. John Wiley and Sons, Inc., 2nd ed., 1992. [4] H. Y. Sutarto, R. A. S., and C. P. P., “Real-time simulator as a test-bed for active gust load alleviation system for an aeroelastic system subject to turbulence.,” (in Indonesian), Internal Report, Vibration Lab., ITB. [5] R. Isermann, J. Schaffnit, and S. Sinsel, “Hardware-in-the-loop simulation for the design and testing of engine-control systems,” Control Engineering Practice, vol. 7, pp. 643—653, 1999. Fig. 2. Components of a simple HIL Simulator

7

Fig. 4. Simulation Modelling Procedure

Fig. 3. Elements of Model Software

Fig. 5. GPS Model

8

Fig. 6. Kalman Filter Algorithm

Fig. 7. Input to trim

Fig. 8.

Fig. 9.

Fig. 10.

Fig. 11.

9

Fig. 12.

Fig. 13.

Fig. 14.