当前位置：首页 > news >正文

ODE45函数——中间变量提取，时变量参数，加速仿真以及运行进度条

news 2025/1/11 0:11:40

System define

The example system:
$G(s)=\frac{U(s)}{s^3-6s^2+11s-6}$
The state-space model:
$\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}$

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);
endplot(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 ¹:

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 2006global odeprogglobvarif 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
elseswitch(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);endend
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 secondsif(t<0)t=-t;end
elseif(length(t1)==1)t=t1;
elset=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'];
elseS=[num2str(t) 's'];
end

The function of odeabort¹:

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];