With the evolution of modern control processes and theories, Model Predictive control has been one of the most widely used control techniques in online applications due to the power of modern computers. MPC (Model Predictive Controller) handles many process control constraints systematically. This control technique is regarded as one of the most advanced techniques for industrial processes.

MPC’s can either be implemented in discrete or continuous time form. However, in this report the MPC has been implemented in the discrete time form as the response time for discrete-time case is a bit faster than the continuous-time case. Hence, the performance is slightly better in the discrete-time case. In this report control strategies for UAV’s (Unmanned Aerial Vehicles) have been discussed.

A UAV or Unmanned Aerial Vehicle is an aircraft with no pilot/crew member/passenger. In recent years, their usage has become more popular in various avenues. Many countries are spending their resources in using these machines in military operations. These vehicles are widely used in surveillance and targeting in war zones. Observations, Agriculture, TV Broadcasting and photography are some of the civil applications of these vehicles. UAV’s are very cost effective when compared with manned surveillance techniques.

The report starts from giving some insights about the history of UAV’s and flight dynamics. In the next part, the equations of motion and arrangement of forces are used in developing a state space model. This is followed by designing a Linear Quadratic Regulator. The same model will then be controlled by designing a MPC.

History and Development of UAV’s

UAV’s have a history touching the 19^th century. In 1848, talented duo of John String Fellow and William Henson built a steam powered propeller driven model UAV. It was known as the ‘Aerial Steam Carriage’ and had a wingspan of 10 foot. In 1896, another researcher named Samuel Langley from the United States developed the Aerodrome-5. In countries like Japan, researchers developed a conventional hydrogen powered UAV and named it as Fugo Balloons. In the 19^th century, UAV’s developed into practical autonomous systems which can be as small as an insect and as big as an airline. Some examples of the 20^th century include Kettering Aerial Torpedo also called the Kettering bug flied with a range of 75 miles. Nevada-Creech air force base used UAV’s to zoom in to attackers to fire missiles. Currently the UAV’s has found the way into the civilian market due to some remarkable developments in the field of electronics and communication. Micro-electronics has seen a significant development over the past couple of years and this has led to small size UAV’s being used for their affordability and simplicity.

Flight Dynamics

In this section, a brief description about the aircraft dynamics is given. Flight dynamics includes all the necessary information about the physical parameters associated with flying objects like aero planes. For clear understanding of the equations of motion governing the aircraft, it is important to understand the coordinate system in airplanes. The three coordinate systems which describe the equations of motion are:

Body-fixed coordinate system
Earth-fixed coordinate system
Atmosphere-fixed coordinate system

The three general reference frames used in airplanes are:

Body Frame
Earth Frame
Stability Frame

Euler Angles (θ, Ø and Ψ) are used to express body-fixed frame in terms of the earth fixed frame.

Body-Fixed Coordinate System

The body-fixed coordinate frame is also called as the ABC frame. As the name suggests this frame is attached to the aircraft body. The x-axis points forward in the aircrafts plane of symmetry as shown in the figure. The y-axis is normal to the aircraft plane of symmetry or the x-axis and pointing to the right side. The z-axis points downwards. The center of gravity is always assumed to coincide with the origin of the body axis system.

Some common terms (see figure 1) associated with the body-fixed frames are:

Roll(Ø)- Roll is referred to the rotation about x-axis.

Pitch(θ)- Pitch is referred to the rotation about y-axis.

Yaw- (Ψ)- Yaw is referred to the rotation about z-axis.

Earth-fixed Coordinate System

The earth-fixed coordinate system is also called as the NED frame. This is because the axis points in the North, East, and Down directions. The NED frame is used to describe the earth reference frame. N represents the position along the northern axis (X_e). E represents the position along the east axis (Y_e). Z represents the local gravity vector (Z_e). In the Fixed-earth coordinate system, the curvature of the earth is often neglected so the NED frame does not rotate. This is also called as the ‘flat-earth approximation’.

Figure 1 – Body-fixed coordinate system representing roll, yaw, and pitch.

Atmosphere-fixed Coordinate System

This axis is defined with respect to the other earth-fixed coordinate system. All the three axes in this system are parallel to the earth-fixed coordinate system. The difference between the two systems is that the atmosphere-fixed coordinate axis moves at a constant velocity with respect to the earth-fixed coordinate axis. This is mainly used in describing the equation of motion.

Apart from the coordinate systems, there are some other terms which are quite often used in describing the dynamics around an airplane. One such term is the Angle of Attack (α) (see figure–). It is the angle between the imaginary line that connects the leading edge and trailing edge of the wing and the relative velocity vector. The angle of attack solely depends on where the aircraft is pointed and to which direction it is moving. It is to be noted that the angle of attack is zero when the stability frame coincides with the body frame.

Figure 2 : Angle of attack (α)

All aerodynamic calculations are based on the body frame. Hence, we need to convert the forces expressed in stability axis to the body axis. The following formulae are used to convert the forces expressed in stability axis to the body axis. Lift and drag forces are converted to the body frame below:

Mw=(w*Rho*S*c)*(Cm0+Cm_alpha*alpha+Cm_delta_e*delta_e)/Iy+Rho*S*c*Cm_alpha*u/(2*Iy);

Mq=(Rho*Va*S*(c^2)*Cm_q)/2*Iy;

Mdelta_e=(Rho*(Va^2)*S*c*Cm_delta_e)/(2*Iy);

%State space representation of longitudinal axis

A_long=[Xu Xw Xq+w -g*cos(theta);Zu Zw Zq-w -g*sin(theta);Mu Mw Mq 0;0 0 1 0];

B_long=[Xdelta_e Xdelta_t;Zdelta_e 0;Mdelta_e 0;0 0];

C_long=eye(4);

D_long=zeros(4,2);

————————————————————–

%Calculating Stability derivatives – Lateral Axis

m=0.568; %Aircraft mass

Q=80.05; % Dynamic Pressure

S=0.2712; %Wing area

Iy=0.11995; %Inertia-x

Iz=0.26547; %Inertia-y

Ix=0.14641; %Inertia-x

C_bar=0.39;%Mean Chord Length

u0=11.432; %(Aircraft Velocity)

g=9.8;

Sw_over_s=0.0187;

e=0.95;

AR=2.42;

tau=0.3;

b=0.81%wing span

theta0=5*pi/180;

%Finding the required Aerodynamic Derivatives:

Cy_beta=-0.06972;

Cy_p=0.0539;

Cy_r=0;

Cy_delta_a=0;

Cn_beta=0.001002;

Cn_p=0.0854;

Cn_r=-0.0154;

Cn_delta_a=-0.0069;

Cl_beta=-0.093;

Cl_p=-0.256;

Cl_r=0.0545;

Cl_delta_a=0.1318;

%Stability deriatives for lateral axis

Y_beta=(Q*S*Cy_beta)/m;

Y_p=(Q*S*b*Cy_p)/(2*u0*m);

Y_r=(Q*S*b*Cy_r)/(2*u0*m);

Y_delta_a=(Q*S*Cy_delta_a)/m;

Y_delta_r=0;

N_beta=(Q*S*b*Cn_beta)/Iz;

N_p=(Q*S*b*b*Cn_p)/(2*u0*Ix);

N_r=(Q*S*b*b*Cn_r)/(2*u0*Ix);

N_delta_a=(Q*S*b*Cn_delta_a)/Iz;

N_delta_r=0;

L_beta=(Q*S*b*Cl_beta)/Ix;

L_p=(Q*S*b*b*Cl_p)/(2*u0*Ix);

L_r=(Q*S*b*b*Cl_r)/(2*u0*Ix);

L_delta_a=(Q*S*b*Cl_delta_a)/Ix;

L_delta_r=0;

%State Space Representation for the Lateral Axis

A_lat=[Y_beta/u0 Y_p -(u0-Y_r) g*cos(theta0);

L_beta/u0 L_p L_r 0;

N_beta/u0 N_p N_r 0;

0 1 0.121 0]

B_lat=[0 Y_delta_r;

L_delta_a L_delta_r;

N_delta_a N_delta_r;

0 0;]

C_lat=eye(4)

D_lat=zeros(4,2)

————————————————————–

%Checking Controllability

Ucont2=[B_lat A_lat*B_lat A_lat^2*B_lat A_lat^3*B_lat]%controllability matrix

unco2 = length(A_lat) – rank(Ucont2) %Checking the number of uncontrollable states

%same code can be used for longitudinal axis by changing the state space matrices.

————————————————————–

%Finding the LQR gain- Longitudinal Axis

%Defining Values of state space equation

C_long_trim=[-0.1193 0.9929 0 -9.7350;

0.9929 0.1193 0 0];

%Trim conditions applied to the output matrix

Q_long=C_long_trim’*C_long_trim

R_long=diag([5 0.1]); %Trim Conditions applied

Skal=B_long*inv(R_long)*B_long’;

%Solving ARE for Pe;

Pe = are(A_long,Skal,Q_long)

%Kalman Gain

Ke=inv(R_long)*B_long’*Pe

————————————————————–

%LQR Gain- Lateral Axis

%Defining Values of state space equation

C_lat_trim=[0 0 0 1];

%Trim conditions applied to the output matrix

Q_lat=C_lat_trim’*C_lat_trim

R_lat=5000; %Trim Conditions applied

Skal1=B_lat*inv(R_lat)*B_lat’;

%Solving ARE for Pe2;

Pe2 = are(A_lat,Skal1,Q_lat)

%Kalman Gain

Ke2=inv(R_lat)*B_lat’*Pe2

%Response for LQR

%Change values for lateral and longitudinal axis

T=0:0.01:10; %Arbitrary sampling instant

U=ones(1001,2);

X0= [9.7 1.2 0 6.9]

sys=ss(A_long-(B_long*Ke),B_long,C_long,D_long);%Change values for lateral and lonigtudinal axis

lsim(sys,U,T,X0)

title(‘response to non zero initial condition with LQR controller’)

————————————————————–

%_Calculating Block toeplitz matrix and then calculating mpc gain values- Longitudinal Axis_

[p_long n_long]=size(C_long);

[n_long,m_long]=size(B_long);

Np_long=10;

Nc_long=4;

v_long=[];

for k=1:Np_long

v_long=[v_long;C_long*A_long^(k-1)*B_long];

end;

F_long=[];

for k=1:Np_long

F_long=[F_long;C_long*A_long^k];

end;

Np_long=10;

Nc_long=4;

Phi_long=zeros(p_long*Np_long, m_long*Nc_long);

%_____________________________________________

% Size of first column of v as well as phi = p X Np

% Size of F = pXNp by n

%__________________________________________________________________

Phi_long=zeros(p_long*Np_long, m_long*Nc_long);

Phi_long(:,1:m_long)=v_long;

for i=2:Nc_long

%Phi(:,i)=[zeros(i-1,1);v(1:Np-i+1,1)];

Phi_long(:,m_long*(i-1)+1:i*m_long)=[zeros(p_long*(i-1),m_long);v_long(1:p_long*(Np_long-i+1),1:m_long)];

end

Rbars_long=F_long(1:p_long*Np_long,n_long) % R is the last coloumn of F

phi_phi_long=Phi_long’*Phi_long

phi_F_long=Phi_long’*F_long

phi_R_long=Phi_long’*Rbars_long

————————————————————–

%_Calculating Block toeplitz matrix and then calculating mpc gain values (Lateral Axis)

[p_lat n_lat]=size(C_lat);

[n_lat,m_lat]=size(B_lat);

Np_lat=10;

Nc_lat=10;

v_lat=[];

for k=1:Np_lat

v_lat=[v_lat;C_lat*A_lat^(k-1)*B_lat];

end;

F_lat=[];

for k=1:Np_lat

F_lat=[F_lat;C_lat*A_lat^k];

end;

Np_lat=10;

Nc_lat=10;

Phi_lat=zeros(p_lat*Np_lat, m_lat*Nc_lat);

%_____________________________________________

% Size of first column of v as well as phi = p X Np

% Size of F = pXNp by n

%__________________________________________________________________

Phi_lat=zeros(p_lat*Np_lat, m_lat*Nc_lat);

Phi_lat(:,1:m_lat)=v_lat;

for i=2:Nc_lat

%Phi(:,i)=[zeros(i-1,1);v(1:Np-i+1,1)];

Phi_lat(:,m_lat*(i-1)+1:i*m_lat)=[zeros(p_lat*(i-1),m_lat);v_lat(1:p_lat*(Np_lat-i+1),1:m_lat)];

end

Rbars_lat=F_lat(1:p_lat*Np_lat,n_lat) % R is the last coloumn of F

phi_phi_lat=Phi_lat’*Phi_lat

phi_F_lat=Phi_lat’*F_lat

phi_R_lat=Phi_lat’*Rbars_lat

————————————————————–

References

Bagheri, S. (2014). Modeling, Simulation and Control System Design for Civil Unmanned Aerial Vehicles. UMEA.
Fahlstrom, P., & Gleason, T. (2012). Introduction to UAV Systems. John Wiley.
Kar, I., & Behera, L. (n.d.). Intelligent Systems and Control. Oxford University Press.
Raemaekers, A. (2007). Design of Model Predictive Controller to Control UAV’s.
Wang, L. (2008). Model Predictive Control System Design and Implentation using MATLAB. Springer.

Index

Air Density 22

Alegbraic Ricatti Equation 27

Angle of Attack 7

Atmosphere-fixed Coordinate System 7

Body-fixed coordinate system 6, 7

Body-Fixed Coordinate System 6

Constraints 33

Control Objectives 33

Controllability 3, 26, 27, 48

cost index 27

Cost or Error Function 34

Elevator Deflection 22

Engine Constant 23

Equations of Motion 3, 9, 10

Euler angles 8

Euler-LaGrange method 9

Flight Dynamics 6

Gravitational constant 23

History 5

Inertia 23, 46

Lauguerre functions 44

Linear Quadratic Regulator (LQR) 3, 27

lsim 29

MIMO 34, 44

Moving Horizon Window 34

Newton’s Approach 9

Pitch Angle 22, 31

Prediction Horizon 34

Receding Horizon Control 34

rotation matrix 8

SISO 44

Slip angle 8

System Dynamics 33

Toeplitz Matrix 35

trim conditions 10, 26, 29, 31

Uncontrollability 26

white noise 35

Wing area 22, 45, 46

Wing Chord 23

Control System Design for Unmanned Aerial Vehicle by Linear Quadratic Regulator and Basic Model Predictive Controller

Control System Design for Unmanned Aerial Vehicle by Linear Quadratic Regulator and Basic Model Predictive Controller

Introduction

History and Development of UAV’s

Flight Dynamics

Body-Fixed Coordinate System

Earth-fixed Coordinate System

Atmosphere-fixed Coordinate System

References

Index

Professor