ODE45函数——中间变量提取，时变量参数，加速仿真以及运行进度条
System define
The example system:
G
(
s
)
=
U
(
s
)
s
3
−
6
s
2
+
11
s
−
6
G(s)=\frac{U(s)}{s^3-6s^2+11s-6}
G(s)=s3−6s2+11s−6U(s)
The state-space model:
d
x
1
d
t
=
x
2
d
x
2
d
t
=
x
3
d
x
3
d
t
=
6
x
1
−
11
x
2
+
6
x
3
+
u
\begin{aligned} \frac{dx_1}{dt} &= x_2 \\ \frac{dx_2}{dt} &= x_3 \\ \frac{dx_3}{dt} &= 6x_1 -11 x_2 + 6x_3 + u \end{aligned}
dtdx1dtdx2dtdx3=x2=x3=6x1−11x2+6x3+u
Section 1: Basic usage (A unit step response)
Note: By limiting the points of time span, the simulation speed can be fastened
tspan = [0,5]; % Simulation time span
% tspan = linspace(0,5,100); % Simulation time span with limited points
y0 = [0;0;0]; % Initial value
[t, y] = ode45(@(t, y) sysOdeFunc(t, y), tspan, y0);
plot(t,y(:,1))
function dydt = sysOdeFunc(t,y)
%--- Unpack variables ---
x1 = y(1);
x2 = y(2);
x3 = y(3);
%--- Controller define ---
u = 1;
%--- Disturbance define ---
d = 1;
%--- System define ---
dx1dt = x2;
dx2dt = x3;
dx3dt = 6*x1 - 11*x2 + 6*x3 + u + d;
%--- Ode array output ----
dydt = [dx1dt;dx2dt;dx3dt];
end
Section 2: Argument dependent on time
Case 1: Implicit form (Only sampled value is known)
tspan = [0,5]; % Simulation time span
% tspan = linspace(0,5,100); % Simulation time span with limited points
y0 = [0;0;0]; % Initial value
ft = linspace(0,5,25); % Here simulate a sampled value and assume it is governed by f(t)
f = ft.^2 - ft - 3; % Here simulate a sampled value and assume it is governed by f(t)
[t, y] = ode45(@(t, y) sysOdeFunc(t, y, ft, f), tspan, y0);
plot(t,y(:,1))
function dydt = sysOdeFunc(t, y, ft, f)
%--- Unpack variables ---
x1 = y(1);
x2 = y(2);
x3 = y(3);
%--- Controller define ---
u = 1;
%--- Disturbance define ---
d = interp1(ft, f, t)
%--- System define ---
dx1dt = x2;
dx2dt = x3;
dx3dt = 6*x1 - 11*x2 + 6*x3 + u + d;
%--- Ode array output ----
dydt = [dx1dt;dx2dt;dx3dt];
end
Case 2: Explicit form (Function expression is known)
tspan = [0,5]; % Simulation time span
% tspan = linspace(0,5,100); % Simulation time span with limited points
y0 = [0;0;0]; % Initial value
d = @(t)sin(t);
[t, y] = ode45(@(t, y) sysOdeFunc(t, y,d), tspan, y0);
plot(t,y(:,1))
function dydt = sysOdeFunc(t,y,d)
%--- Unpack variables ---
x1 = y(1);
x2 = y(2);
x3 = y(3);
%--- Controller define ---
u = 1;
%--- Disturbance define ---
% No need
%--- System define ---
dx1dt = x2;
dx2dt = x3;
dx3dt = 6*x1 - 11*x2 + 6*x3 + u + d(t);
%--- Ode array output ----
dydt = [dx1dt;dx2dt;dx3dt];
end
Section 3: Output additional variable
tspan = [0,5]; % Simulation time span
% tspan = linspace(0,5,100); % Simulation time span with limited points
y0 = [0;0;0]; % Initial value
d = @(t)sin(t);
[t, y] = ode45(@(t, y) sysOdeFunc(t, y, d), tspan, y0);
%--- Additional variables extraction ----
for k = 1:max(size(t))
[~, uout(k,:)] = sysOdeFunc(t(k), y(k,:), d);
end
plot(t,y(:,1)), title('System state')
figure
plot(t,uout), title('Additional variable')
function [dydt,u] = sysOdeFunc(t,y,d)
%--- Unpack variables ---
x1 = y(1);
x2 = y(2);
x3 = y(3);
%--- Controller define ---
u = sin(t);
%--- Disturbance define ---
% No need
%--- System define ---
dx1dt = x2;
dx2dt = x3;
dx3dt = 6*x1 - 11*x2 + 6*x3 + u + d(t);
%--- Ode array output ----
dydt = [dx1dt;dx2dt;dx3dt];
end
Section 4: Add the ODE process and abort selection
The layout would be like below:
The code is as below
tspan = [0,5]; % Simulation time span
% tspan = linspace(0,5,100); % Simulation time span with limited points
y0 = [0;0;0]; % Initial value
odeopts = odeset('OutputFcn', @odeprog, 'Events',@odeabort);
[t, y] = ode45(@(t, y) sysOdeFunc(t, y), tspan, y0, odeopts);
plot(t,y(:,1))
function dydt = sysOdeFunc(t,y)
%--- Unpack variables ---
x1 = y(1);
x2 = y(2);
x3 = y(3);
%--- Controller define ---
u = 1;
%--- Disturbance define ---
d = 1;
%--- System define ---
dx1dt = x2;
dx2dt = x3;
dx3dt = 6*x1 - 11*x2 + 6*x3 + u + d;
%--- Ode array output ----
dydt = [dx1dt;dx2dt;dx3dt];
end
Two functions are defined as follows:
The function of odeprog 1:
function status = odeprog(t,y,flag,varargin)
%status = odebarplot(t,y,flag,varargin)
% ODE progress display function with interrupt control
% Displays a vertical bar plot that fills as the simulation
% nears its completion. Also displays time ellapsed and estimates
% time remaining in the simulation. To avoid computation burden
% refreshes are limited to every 0.5 seconds.
%
% Tim Franklin
% Virginia Tech
% Jesse Norris
% Wake Forrest
% May 2006
global odeprogglobvar
if nargin < 3 || isempty(flag)
if(etime(clock,odeprogglobvar(8:13))>0.5)
tfin=odeprogglobvar(1);
sstrt=odeprogglobvar(2:7);
figure(95);
perc=t(end)/tfin;
area([t(end) tfin-t(end);t(end) tfin-t(end)]);
title([num2str(perc*100) '%']);
set(findobj('Tag','eltime'),'String',etimev(clock,sstrt));
set(findobj('Tag','esttime'),'String',etimev(etime(clock,sstrt)/perc*(1-perc)));
odeprogglobvar(8:13)=clock;
end
else
switch(flag)
case 'init'
odeprogglobvar=zeros(1,13);
odeprogglobvar(1)=t(end);
odeprogglobvar(2:7)=clock;
odeprogglobvar(8:13)=clock;
tfin=odeprogglobvar(1);
sstrt=odeprogglobvar(2:7);
figure(95);
set(gcf,'Position',[4,40,100,500]);
axes('Position',[0.5,0.25,0.25,0.6]);
axis([1,2,0,tfin]);
set(gca,'XTickLabel',[],'NextPlot','replacechildren');
ylabel('Simulation Progress - Time (s)');
title('0%');
area([0 tfin;0 tfin]);
uicontrol('Style', 'pushbutton', 'String', 'Abort','Position', [7 460 90 30], 'Callback', 'close(gcf)')
uicontrol('Style', 'text', 'String', 'Ellapsed Time','Position', [7 100 90 15])
uicontrol('Style', 'text', 'Tag', 'eltime', 'String', etimev(clock,sstrt),'Position', [7 80 90 15])
uicontrol('Style', 'text', 'String', 'Time Remaining','Position', [7 60 90 15])
uicontrol('Style', 'text', 'Tag', 'esttime', 'String', num2str(inf),'Position', [7 40 90 15])
pause(0.1);
case 'done'
if(ishandle(95))
close(95);
end
end
end
status = 0;
drawnow;
function [S] = etimev(t1,t0)
%ETIMEV Verbose Elapsed time.
% ETIMEV(T1,T0) returns string of the days, hours, minutes, seconds that have elapsed
% between vectors T1 and T0. The two vectors must be six elements long, in
% the format returned by CLOCK:
%
% T = [Year Month Day Hour Minute Second]
% OR
% ETIMEV(t), t in seconds
if(exist('t1')&exist('t0')&length(t1)>2&length(t0)>2)
t=etime(t1,t0); %Time in seconds
if(t<0)
t=-t;
end
elseif(length(t1)==1)
t=t1;
else
t=0;
end
days=floor(t/(24*60*60));
t=t-days*24*60*60;
hours=floor(t/(60*60));
t=t-hours*60*60;
mins=floor(t/60);
t=floor(t-mins*60);
if(days>0)
S=[num2str(days) 'd ' num2str(hours) 'h ' num2str(mins) 'm ' num2str(t) 's'];
elseif(hours>0)
S=[num2str(hours) 'h ' num2str(mins) 'm ' num2str(t) 's'];
elseif(mins>0)
S=[num2str(mins) 'm ' num2str(t) 's'];
else
S=[num2str(t) 's'];
end
The function of odeabort1:
function [value,isterminal,direction]=odeabort(t,S,varargin)
%Other Events Set Here...ie:
% value(2)=max(abs(S(1:18)))-pi/2;
%Test to see if 'simstop' box is closed
value(1)=double(ishandle(95));
isterminal=[1];
direction=[0];
https://www.mathworks.com/matlabcentral/fileexchange/9904-ode-progress-bar-and-interrupt ↩︎ ↩︎