[Quadrotor] Simulate the translational and rotational dynamics of a quadrotor (Simulink simulation implementation)

Welcome to this blog

Advantages of bloggers:Blog content should be as thoughtful and logical as possible for the convenience of readers.

Motto:He who travels a hundred miles is half as good as ninety.

The directory of this article is as follows:

Table of Contents

1 Overview

2 Operation results

3 References

4 Matlab code, Simulink simulation and literature


1 Overview

Abstract-Due to the growing interest in unmanned aerial vehicles (UAVs) for both civilian and military applications, research on autonomous micro-flying robots has been significantly enhanced. This article summarizes the final results of the modeling and control portion of the OS4 project, which focused on quadcopter design and control. The article introduces a simulation model that takes into account changes in aerodynamic coefficients caused by aircraft motion. The control parameters found using this model were successfully applied to the helicopter without readjustment. The final part of this paper describes the control method (integral feedback) and our proposed complete control scheme for the quadcopter (attitude, altitude and position). Finally, the results of autonomous takeoff, hovering, landing and obstacle avoidance are presented.

Flying objects have always had a great attraction for mankind, encouraging various research and development. The project began in 2003, when the robotics community showed increasing interest in the development of Unmanned Aerial Vehicles (UAVs). The scientific challenges of designing and controlling drones in crowded environments and the lack of existing solutions are highly motivating. On the other hand, the wide range of application areas in the military and civilian markets encourages financial support for UAV-related projects. From the beginning of the project it was decided to focus on one specific configuration: the quadcopter. This choice stems not only from its dynamic properties, which represent an attractive control problem, but also from design issues. Integrating sensors, actuators and intelligence into a lightweight vertical flight system while maintaining good uptime is no simple matter.

A. Current situation
The current state of quadcopter control has changed dramatically over the past few years. The number of projects addressing this problem has increased significantly. Most of these projects are based on commercial toys, such as the Draganflyer [1], which were later modified to add more sensing and communication capabilities. Only a few teams have solved the MFR design problem. The paper [2] lists some of the most important quadcopter projects in the past 10 years. The Mesicopter project [3] started in 1999 and ended in 2001. It aims to study the feasibility of centimeter-scale quadcopters. E. Altug introduced dual-camera visual feedback control in 2003 in his paper [4].

See Part 4 for detailed explanation.

2 Operation results

Matlab code:

% KF_setup.m 19/11/2013
%Quadrotor Sim
%
%
% Purpose: to declare initial values
%
%%
%% sensor noises
position_uncertainty_var = (20/3600*pi/180)^2*ones(3,1);
%% simulation set up
step_time = 0.5; % simulation step time(sec)
end_time = 1000; % simulation end time (sec)
%end_time = 86400;
%% attitude estimator gains
Tatd = 0.5; % attitude estimator update time (sec)
Tqint = 0.5; % discrete quaternion integration period (sec)
Tsen_out = 0.5; % sensor output period (sec)
TkfProp = 0.5; % Kalman filter propagation period (sec)
KfupdatePeriodInCycle = 1; % Kalman filter update period (propagation cycle)
f_bw_atd = 0.02; % attitude determination bandwidth (hz)
%f_bw_atd = 0.005;
zeta = 0.7;
Krp = (2*pi*f_bw_atd)^2 * eye(3);
Kpp = 2*zeta*2*pi*f_bw_atd*eye(3);
qest0 = [0*1e-4; 0; 0; 1]; % initial estimator quaternion
delta_west0 = zeros(3,1); % initial deviation of estimator angular rate (rad/sec)
max_delta_w = 0.1*pi/180;
delta_w_lim = 2e-4; %0.1/pi/Tqint;
delta_th_lim= 1e-4; %0.1*pi/180/Tqint;
q0 = [0; 0; 0; 1];

%% for estimate error standard deviation prediction calculation
wn=sqrt(diag(Krp));
k=sqrt((wn.^4 + 4*zeta^2)./(4*zeta*wn));
%% for using Lyapunove equation to solve for expected estimation error
C=[1 0]; K=[Kpp(1,1);Krp(1,1)]; A=[0 1;0 0]-K*C; B=K;
H=[1 0]; K=[Kpp(1,1);Krp(1,1)]*Tatd; F=[1 Tatd;0 1]-K*H; G=K;
%% Kalman filter setups
Fmat = [eye(3) TkfProp*eye(3);zeros(3,3) eye(3)];
Hmat = [eye(3) zeros(3,3)];
therr0 = max([abs(qest0(1:3)); 5*1e-4]); % initial error estimate, assuming q0=[0 0 0 1]
P0 = diag([therr0^2*ones(1,3) 3e-6^2*ones(1,3)]);
R = TkfProp*KfupdatePeriodInCycle*diag(position_uncertainty_var); -3^2*eye(3)*
Q = diag([1e-5^2*ones(1,3), 1e-7^2*ones(1,3)])*TkfProp;
max_rate = pi/180;
P0 = diag([1e-32*ones(1,3) 1e-5^2*ones(1,3)]);
Q = diag([1e-5^2*ones(1,3), 5e-6^2*ones(1,3)])*TkfProp;
max_bias = 1*pi/180/3600;
%% start simulation
Tcapt = Tsen_out; % sim variable capture rate (sec)

% KF_setup.m 19/11/2013
%Quadrotor Sim
%
%
% Purpose: to declare initial values
%
%%
%% sensor noises
position_uncertainty_var = (20/3600*pi/180)^2*ones(3,1);
%% simulation set up
step_time = 0.5; % simulation step time(sec)
end_time = 1000; % simulation end time (sec)
%end_time = 86400;
%% attitude estimator gains
Tatd = 0.5; % attitude estimator update time (sec)
Tqint = 0.5; % discrete quaternion integration period (sec)
Tsen_out = 0.5; % sensor output period (sec)
TkfProp = 0.5; % Kalman filter propagation period (sec)
KfupdatePeriodInCycle = 1; % Kalman filter update period (propagation cycle)
f_bw_atd = 0.02; % attitude determination bandwidth (hz)
%f_bw_atd = 0.005;
zeta = 0.7;
Krp = (2*pi*f_bw_atd)^2 * eye(3);
Kpp = 2*zeta*2*pi*f_bw_atd*eye(3);
qest0 = [0*1e-4; 0; 0; 1]; % initial estimator quaternion
delta_west0 = zeros(3,1); % initial deviation of estimator angular rate (rad/sec)
max_delta_w = 0.1*pi/180;
delta_w_lim = 2e-4; %0.1/pi/Tqint;
delta_th_lim= 1e-4; %0.1*pi/180/Tqint;
q0 = [0; 0; 0; 1];

%% for estimate error standard deviation prediction calculation
wn=sqrt(diag(Krp));
k=sqrt((wn.^4 + 4*zeta^2)./(4*zeta*wn));
%% for using Lyapunove equation to solve for expected estimation error
C=[1 0]; K=[Kpp(1,1);Krp(1,1)]; A=[0 1;0 0]-K*C; B=K;
H=[1 0]; K=[Kpp(1,1);Krp(1,1)]*Tatd; F=[1 Tatd;0 1]-K*H; G=K;
%% Kalman filter setups
Fmat = [eye(3) TkfProp*eye(3);zeros(3,3) eye(3)];
Hmat = [eye(3) zeros(3,3)];
therr0 = max([abs(qest0(1:3)); 5*1e-4]); % initial error estimate, assuming q0=[0 0 0 1]
P0 = diag([therr0^2*ones(1,3) 3e-6^2*ones(1,3)]);
R = TkfProp*KfupdatePeriodInCycle*diag(position_uncertainty_var); -3^2*eye(3)*
Q = diag([1e-5^2*ones(1,3), 1e-7^2*ones(1,3)])*TkfProp;
max_rate = pi/180;
P0 = diag([1e-32*ones(1,3) 1e-5^2*ones(1,3)]);
Q = diag([1e-5^2*ones(1,3), 5e-6^2*ones(1,3)])*TkfProp;
max_bias = 1*pi/180/3600;
%% start simulation
Tcapt = Tsen_out; % sim variable capture rate (sec)

3 References

Some of the content in the article is quoted from the Internet. The source will be indicated or cited as a reference. It is inevitable that there will be some unfinished information. If there is anything inappropriate, please feel free to contact us to delete it.

[1] Yang Wei, Zhou Bingqian, Wu Zhigang, et al. Design and simulation analysis of quad-rotor aircraft [J]. Computer CD-ROM Software and Applications, 2014, 17(16):2.DOI:CNKI:SUN:GPRJ.0.2014-16 -154.

[2] Gao Yan, Yu Dan. Modeling and control algorithm simulation of quad-rotor aircraft [J]. Industrial Control Computer, 2014(9):3.DOI:10.3969/j.issn.1001-182X.2014.09.045.

[3] Qiao Weiwei. Research and simulation of quadcopter flight control system [D]. North China University, 2012. DOI: 10.7666/d.D316360.