solution of ode's with inputs using ode45

14 views (last 30 days)
savvas
savvas on 23 Oct 2015
Commented: savvas on 23 Oct 2015
Why do these two pieces of script give the same output?
Linear Dynamics:
% This code calculates the solution to the Linear system
%-------------------------------------------------------%
function LinSys = LinSol()
tspan = [0 10]; % time span
iniCon = [-1; -1; 0; 0; 0; 0]; % initial conditions
[t,y] = ode45(@(t,x) sys(t,x),tspan,iniCon);
plot(t,y,'linewidth',2);
legend({'\eta_1';'\eta_2';'\eta_3';'\eta_4';'\eta_5';'\eta_6'})
title('Linear Response');
xlabel('Time (s)')
ylabel('Position (m)');
grid on;
% This function simulates the Linear System
%-----------------------------------------%
function dx = sys(t,x)
m = 4; % mass
J = 0.0475; % moment of inertia
r = 0.25; % distance from c.o.m to vector force
g = 9.8; % gravity
c = 0.05; % damping coefficient
A = [0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
0 0 -g -c/m 0 0
0 0 0 0 -c/m 0
0 0 0 0 0 0];
B = [0 0
0 0
0 0
1/m 0
0 1/m
r/J 0];
K = findLQR;
u = -K*x; % control law
% Linear System
%-------------%
dx = A*x + B*u;
% This function calculates the gain matrix K
%------------------------------------------%
function K = findLQR()
p = 0.1;
Q=diag([1 1 1 1 1 1]); % state cost matrix
R=diag([p p]); % input cost matrix
[K,~,~]=lqr(A,B,Q,R); % optimal gain matrix
end
end
end
and Nonlinear Dynamics:
% This code calculates the solution to the Nonlinear system
%---------------------------------------------------------%
function LinSys = NonLinSol()
clearvars
tspan = [0:0.01:10]; % time span
iniCon = [-1; -1; 0; 0; 0; 0]; % initial conditions
[t,y] = ode45(@(t,x) sys(t,x),tspan,iniCon);
plot(t,y,'linewidth',2);
legend({'\eta_1';'\eta_2';'\eta_3';'\eta_4';'\eta_5';'\eta_6'})
title('Nonlinear Response');
xlabel('Time (s)')
ylabel('Position (m)');
grid on;
% This function simulates the Nonlinear System
%--------------------------------------------%
function dx = sys(t,x)
m = 4; % mass
J = 0.0475; % moment of inertia
r = 0.25; % distance from c.o.m to vector force
g = 9.8; % gravity
c = 0.05; % damping coefficient
A = [0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
0 0 -g -c/m 0 0
0 0 0 0 -c/m 0
0 0 0 0 0 0];
B = [0 0
0 0
0 0
1/m 0
0 1/m
r/J 0];
K = findLQR;
dx = zeros(6,1);
u = -K*x; % control law
% Nonlinear Dynamics
%------------------%
dx(1) = x(4);
dx(2) = x(5);
dx(3) = x(6);
dx(4) = -g*sin(x(3))-(c/m)*x(4)+(u(1)/m)*cos(x(3))-(u(2)/m)*sin(x(3));
dx(5) = g*cos(x(3))-g-(c/m)*x(5)+(u(1)/m)*sin(x(3))+(u(2)/m)*cos(x(3));
dx(6) = (r/J)*u(1);
% This function calculates the gain matrix K
%------------------------------------------%
function K = findLQR()
p = 0.1;
Q=diag([1 1 1 1 1 1]); % state cost matrix
R=diag([p p]); % input cost matrix
[K,~,~]=lqr(A,B,Q,R); % optimal gain matrix
end
end
end

Answers (1)

@Johannes
@Johannes on 23 Oct 2015
Edited: @Johannes on 23 Oct 2015
Hello,
i agree with your linear approach. But in your nonlinear script you are mixing linear control theory with a nonlinear system. You calculate the matrix K by the matrices A and B of the linear System. But you have to use a nonlinear approach to control the nonlinear system.
I have modified your code and there you will see that your results are not equal.
Best Regards, Johannes
  3 Comments
@Johannes
@Johannes on 23 Oct 2015
You suggested that both scripts give the same output. I just wanted to show you that this is not true. If you zoom into the plot you will see that there are differences between the lines. The x x x lines represent the linear system.
There are some nice files from the ETHZ about Nonlinear control. http://control.ee.ethz.ch/~apnoco/lectures.php
Unfortunately i can not recommend you an approach. This depends on what you are trying to achieve.
Best Regards,
Johannes
savvas
savvas on 23 Oct 2015
Yes I know that there is a small difference but still the response is not the expected one for the nonlinear system. I have the plot of the nonlinear system using a much different approach to get it and it is the correct one. Thanks for your input.

Sign in to comment.

Categories

Find more on Systems of Nonlinear Equations in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!