MATLAB/Simulink利用状态空间方程求解动力学方程组的问题求助

[weixin_46859090]

1.首先，我需要求解一个具有14个广义坐标的动力学方程组，初步建立Simulink模型，利用状态空间方程求解，如下所示，求解成功；

 

2.接下来我需要对该模型进行调频控制，因此暂时利用傅里叶变换求解其主结构的频率，由所得的频率来改变减振装置的刚度，根据刚度值计算得到控制力进行计算，求解成功，模型如下所示：

不过与步骤1中得到的数据对比发现，该调频控制几乎没有任何效果，与被动控制的动力学响应几乎没有差别。

思考过后，发现是状态空间方程的问题，据上网了解，State space模块输入的矩阵一般读取的是由脚本文件中输入到工作区的矩阵，且无法实时改变。但我要做的工作是根据得到的频率实时更新刚度值，状态空间方程中所涉及的刚度矩阵也需要实时变化，因此考虑不使用State space模块，而是自己利用MATLAB Function模块自己编写方程进行求解。

3.对步骤2中的模型进行小修，只是将State space模块替换成自己编写的MATLAB Function模块，如下所示。

MATLAB Function模块，即上图的State space模块的代码如下所示：

 

function dis = state_space(t, u, k_new, MGain, KGain, CGain)
% 输入:
%   t: 时间向量
%   u: 外载荷矩阵 (每列对应一个时间步的载荷向量)
%   k_new: 新的TMD刚度值 (标量)
% 输出:
%   dis: 位移响应矩阵 (每行对应一个自由度，每列对应一个时间步)

M=MGain;
C=CGain;
K0=KGain;

% 获取系统参数
n_dof = size(M, 1); % 自由度数量
n_steps = length(t); % 时间步数

% 初始化输出矩阵
dis = zeros(n_dof, n_steps);
velocity_response = zeros(n_dof, n_steps);

% 更新刚度矩阵中的TMD部分
% 假设TMD刚度位于K矩阵的最后两个对角线元素
K = K0;
K(end-1:end, end-1:end) = [k_new, 0; 0, k_new];

% 初始条件 (假设初始位移和速度均为零)
x0 = zeros(2*n_dof, 1);

% 计算状态空间矩阵
A = [zeros(n_dof), eye(n_dof); -M\K, -M\C];
% B = [zeros(n_dof, n_dof*2); M\eye(n_dof), eye(n_dof)];
B = [zeros(n_dof, size(u, 1)); M\eye(n_dof)];

% 使用数值积分方法求解状态方程
% 这里使用简单的欧拉法，您可以根据需要更改为更精确的方法
dt = 0.01; % 时间步长
x = x0;

for i = 1:n_steps
    % 计算状态导数
    dx = A * x + B * u;

    % 更新状态
    x = x + dx * dt;

    % 存储响应
    dis(:, i) = x(1:n_dof);
    velocity_response(:, i) = x(n_dof+1:end);
end

end

求解完成，没有报错，但出现了两个问题：

（1）最关键的问题，我求解的Dis.文件，即位移响应全变成了0，查看整个大矩阵，除了时间序列，剩余14个广义坐标的位于响应序列均为0；

（2）最初步骤1时我的外载荷包含地震载荷的，后续暂时不需要地震载荷就注释掉了，但mux模块也注释的话会发生看不懂的报错，因此作为半吊子的我本着代码能动就不动他的原则未作处理；而在步骤3时，不断尝试中发现mux模块需要注释掉，否则也会报错，simulink中的报错信息如下，不知道是什么原因，也不知道是否是导致问题（1）的罪魁祸首。

要串联的数组的维度不一致。
函数 'State Space1' (#902.650.690)，行 33，列 5:
"[zeros(n_dof, size(u, 1)); M\eye(n_dof)]"

因此小白请教各位大佬，我在步骤2到步骤3期间，是什么导致求解不出位移响应？该如何修改才能实现这个调频控制呢？