大规模非线性动力系统的高效可识别性、可控性和可观测性检测研究（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客 

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥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

🌈4 Matlab代码实现