👨🎓个人主页:研学社的博客
💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。
⛳️座右铭:行百里者,半于九十。
📋📋📋本文目录如下:🎁🎁🎁
目录
💥1 概述
本文对大规模非线性动力系统的高效可识别性、可控性和可观测性检测研究,并用Matlab代码实现之。详细文章讲解见参考文献。
📚2 运行结果
部分结果图
部分代码:
function [txxp,xxpt] = pxxssen(t,xt,uss,pe,fxup,dt,x0,ix0e,ad,xpi,xxi)
% pxxssen : Steady state sensitivity matrix for identifiability
%
% [txxp,xxpt] = pxxssen(t,xt,uss,pe,fxup,dt,x0,ix0e,ad)
%
% Inputs
% t : time t
% xt : system state response
% uss : steady input
% pe : estimated parameters
% fxup : function handle to system dynamics with pe as argument
% dt : time-step scaling and squaring
% x0 : initial state
% ix0e : indices to be estimated initial states <0 all 0 none
% ad : (optional) if nonzero automatic differentiation else finite differences
% default: automatic differentiation
% xpi : (optional) Initial values state sensitivity to parameters
% default zeros(nx,np)
% xxi : (optional) Initial values state sensitivity to parameters
% default eye(nx,nx)
%
% Outputs
% txxp : time axis state sensitivity matrix xxpt
% xxpt : state sensitivity matrix xxpt
%
% GvW 14-11-2018
%% Check inputs
if nargin<6; error(' To few inputs');
elseif nargin<8; ix0e=[]; x0=[]; ad=1;
elseif nargin<9; ad=1;
elseif nargin<10; xpi=[];
elseif nargin<11; xxi=[];
end % Default: automatic differentiation
th=t(1); xh=xt(1,:)'; % Extract initial time and system state
if size(uss,2)>1; uss=interp1(uss(:,1),uss(:,2:end),th)'; end
np=length(pe); % Number of estimated parameters
[nt,nx]=size(xt); % State dimension nx
if ~isempty(x0) && nx~=length(x0); error(' xt,x0 incompatible'); end
nu=length(uss);
%% Initialize sensitivity matrix
if isempty(xxi);
xxi=zeros(nx,length(ix0e)); for ih=ix0e; xxi(ih,ih)=1; end
end
if max(ix0e)<1;
xph=xpi; % Initial state sensitivity if x0 not estimated
else
xph=[xpi,xxi]; % Initial state sensitivity if x0 estimated
pe=[pe; x0(ix0e)];
end
txxp=th; %Initialize txxp
xxpt=[xh; xph(:)]'; % Initialize state sensitivity storage
% Empty pe if only x0 is identified
if np; peh=pe; else peh=[]; end;
%% Single df/dx, df/dp computation (because of steady state)
if ~nu
[dfdx]=vfdiff(@(t,xh)fxup(t,xh,peh),th,xh,ad);
if np; [dfdp]=vfdiff(@(t,peh)fxup(t,xh,peh),th,peh,ad);
else dfdp=[]; end
else
[dfdx]=vfdiff(@(t,xh)fxup(t,xh,uss,peh),th,xh,ad);
if np; [dfdp]=vfdiff(@(t,peh)fxup(t,xh,uss,peh),th,peh,ad);
else dfdp=[]; end
end
%% Computation phi(0,dt), gam(0,dt)
ns=ceil(abs((t(end)-t(1))/dt));
dt=(t(end)-t(1))/ns;
if dt>0; A=dfdx; else A=-dfdx'; end
if isempty(dfdp)
phi=expm(A*dt);
else
[phi,gam]=c2d(A,dfdp,dt);
end
%% Integration sensitivity matrix; xh unchanged steady state
for i1=1:ns
xph=phi*xph; if np; xph=xph+gam; end;
th=t(1)+i1*dt; txxp=[txxp;th]; xxpt=[xxpt; [xh;xph(:)]'];
end
🎉3 参考文献
部分理论来源于网络,如有侵权请联系删除。
[1]L.G. Van Willigenburg, J.D. Stiger, J. Molenaar, 2022, "Sensitivity matrices as keys to local structural system properties of large-scale nonlinear system