
简介本资源是一份面向控制工程、智能系统建模与物理信息机器学习方向的MATLAB实践案例专为具备基础动力学知识和MATLAB编程能力的高年级本科生、研究生及科研初学者设计用于掌握如何用物理信息神经网络PINNs求解经典二阶振动系统——质量-弹簧-阻尼器模型。资源共7个文件含3个核心MATLAB脚本用于数据生成、PINNs构建与预测可视化、2张结果对比图png、1个预置训练数据集mat及1份说明文档md包体仅46KB轻量易读、结构清晰便于快速复现与二次开发。目前已有92人学习下载适合希望深入理解PINNs在连续动力系统建模中应用逻辑的学习者。读者可直接运行代码完成从微分方程建模、神经网络架构设计、数据驱动训练到响应预测与误差分析的全流程附带可视化图表与模块化函数显著降低物理神经网络入门门槛。1. 用物理信息神经网络PINN在 MATLAB 中求解质量-弹簧-阻尼器系统不是替代而是增强传统建模能力你手头有一组实测的位移-时间数据系统结构明确单自由度质量块、线性弹簧、粘性阻尼器但参数未知或存在非线性退化或者你想在不依赖大量标注数据的前提下让神经网络输出严格满足牛顿第二定律 $ m\ddot{x} c\dot{x} kx f(t) $ 的解。这时物理信息神经网络PINN不是要取代 Simulink 或ode45而是把先验物理规律“编译”进网络训练过程——它让模型既拟合观测又服从方程约束。MATLAB 用户常误以为 PINN 必须用 PythonPyTorch 实现其实从 R2021b 起Deep Learning Toolbox 已原生支持自定义损失函数、符号微分Symbolic Math Toolbox 配合jacobian/diff和自动微分dlgradient完全可在纯 MATLAB 环境中构建可微分、可验证、可部署的 PINN 求解器。本文面向有 ODE 建模经验、熟悉trainNetwork流程但尚未将物理方程嵌入神经网络的工程师提供一套可直接运行、参数可调、误差可量化、结果可导出的完整实现路径。2. 构建 PINN 求解器从物理方程到可微分网络架构2.1 明确控制方程与边界/初始条件的数学表达质量-弹簧-阻尼器系统的动力学本质是二阶常微分方程ODE $$ m \frac{d^2 x}{dt^2} c \frac{dx}{dt} k x f(t) $$ 其中 $ x(t) $ 是位移响应$ m, c, k $ 是待识别参数或已知常数$ f(t) $ 是外激励如阶跃、正弦或实测力信号。PINN 的核心思想是不直接求解该 ODE而是构造一个神经网络 $ \hat{x}_\theta(t) $使其输出在训练域内同时最小化两类残差数据残差在已知测量点 $ t_i $ 上$ \hat{x}_\theta(t_i) $ 与实测位移 $ x_i $ 的均方误差物理残差将 $ \hat{x}\theta(t) $ 代入 ODE 左侧计算残差函数 $ \mathcal{R}\theta(t) m \ddot{\hat{x}}\theta c \dot{\hat{x}}\theta k \hat{x}\theta - f(t) $并在采样点上最小化 $ |\mathcal{R}\theta|^2 $。提示物理残差必须显式计算二阶导数。MATLAB 中不能仅靠gradient多次近似精度不足且不可导必须使用符号微分或自动微分。本方案采用dlgradientdlarray实现高阶导数兼容 GPU 加速且梯度稳定。2.2 设计可微分网络结构与输入输出接口我们选用深度为 4、每层 50 个神经元的全连接网络dlnetwork激活函数为tanh优于 ReLU因二阶导数连续利于物理残差收敛。输入为标量时间 $ t $输出为标量位移 $ \hat{x}(t) $。关键在于网络输出必须是dlarray类型且所有微分操作在dlarray上进行。% 定义网络层MATLAB R2022b layers [ featureInputLayer(1,Normalization,none) fullyConnectedLayer(50) tanhLayer fullyConnectedLayer(50) tanhLayer fullyConnectedLayer(50) tanhLayer fullyConnectedLayer(1) regressionLayer]; % 构建 dlnetwork 对象启用自动微分 net dlnetwork(layers,OutputNames,{output}); % 初始化权重避免对称性导致训练停滞 net initialize(net);网络本身不包含物理方程物理约束通过后续的自定义损失函数注入。输入 $ t $ 需归一化到 [-1,1] 区间提升训练稳定性例如t_norm 2*(t - t_min)/(t_max - t_min) - 1;。2.3 实现物理残差计算用dlgradient获取高阶导数MATLAB 的dlgradient只支持一阶导数但可通过嵌套调用获得二阶导数。以下函数computePhysicsResidual接收归一化时间t_dldlarray、网络对象net、参数m,c,k和激励f_t返回物理残差 $ \mathcal{R}(t) $function R computePhysicsResidual(t_dl, net, m, c, k, f_t) % 前向传播x_hat net(t) x_hat forward(net, t_dl); % 一阶导数v_hat dx/dt v_hat dlgradient(sum(x_hat), t_dl, RetainData, true); % 二阶导数a_hat d²x/dt²注意t_dl 是归一化时间需链式法则修正 % dt_norm/dt 2/(t_max-t_min)故 d²x/dt² d²x/dt_norm² * (dt_norm/dt)² dt_norm_dt 2/(t_max - t_min); % 需在调用前定义 t_max, t_min a_hat dlgradient(sum(v_hat), t_dl, RetainData, true) * dt_norm_dt^2; % 物理残差m*a c*v k*x - f(t) R m * a_hat c * v_hat k * x_hat - f_t; end注意dt_norm_dt是归一化尺度因子必须在训练循环外预先计算并传入否则dlgradient无法对其求导。若忽略此因子残差量纲错误训练必然发散。f_t应为与t_dl同尺寸的dlarray支持符号表达式如sin(2*pi*t_dl)或插值数组。2.4 组装 PINN 训练循环混合损失与自适应权重总损失函数为加权和 $$ \mathcal{L} \lambda_{\text{data}} \cdot \mathcal{L}{\text{data}} \lambda{\text{physics}} \cdot \mathcal{L}{\text{physics}} $$ 其中 $ \mathcal{L}{\text{data}} \frac{1}{N_d}\sum_i (\hat{x}\theta(t_i) - x_i)^2 $$ \mathcal{L}{\text{physics}} \frac{1}{N_p}\sum_j \mathcal{R}_\theta^2(t_j) $。权重 $ \lambda $ 决定数据拟合与物理守恒的优先级。实践中固定权重易导致某一项主导训练推荐使用自适应策略初期侧重物理残差强制网络学习方程后期提升数据权重精调拟合。以下为带权重调度的训练主干% 初始化优化器Adam opt adamOptimizer(LearnRate, 0.001); % 预分配训练数据t_data, x_data 为实测点t_physics 为物理残差采样点 t_data_dl dlarray(t_data_norm, CB); % C: channel, B: batch x_data_dl dlarray(x_data, CB); t_physics_dl dlarray(linspace(-1,1,200), CB); % 200个物理点 f_physics sin(2*pi*t_physics_dl*0.5); % 示例激励 for epoch 1:1000 % 动态权重物理权重随 epoch 递减数据权重递增 lambda_p max(1.0, 10.0 - epoch/100); lambda_d min(1.0, epoch/500); % 计算数据损失 x_pred predict(net, t_data_dl); loss_data mean((x_pred - x_data_dl).^2); % 计算物理损失调用前述函数 R_physics computePhysicsResidual(t_physics_dl, net, m_true, c_true, k_true, f_physics); loss_physics mean(R_physics.^2); % 总损失 loss_total lambda_d * loss_data lambda_p * loss_physics; % 反向传播更新网络 [gradients, state] dlgradient(loss_total, net.Learnables); net update(net, gradients, opt); opt update(opt, gradients, net.Learnables); % 每100轮打印损失 if mod(epoch,100)0 fprintf(Epoch %d: L_data%.2e, L_physics%.2e\n, ... epoch, double(loss_data), double(loss_physics)); end end提示t_data_norm和t_physics_dl必须同为dlarray且维度匹配CB格式。predict函数自动处理dlarray输入无需手动转换。state用于保存优化器内部状态如 Adam 的动量不可省略。3. 参数识别与系统辨识从 PINN 输出反推未知物理量3.1 将未知参数作为网络可学习变量嵌入 PINN当 $ m, c, k $ 未知时不能将其设为常数传入computePhysicsResidual。正确做法是将参数声明为dlarray并加入网络Learnables使它们与网络权重一同被优化。修改网络初始化% 初始化参数为 dlarray对数空间初始化保证正定 m_init dlarray(log(1.0), U); % U: unformatted scalar c_init dlarray(log(0.5), U); k_init dlarray(log(10.0), U); % 将参数加入 Learnables 列表 net.Learnables [net.Learnables; ... struct(Parameter,m,Value,m_init); ... struct(Parameter,c,Value,c_init); ... struct(Parameter,k,Value,k_init)];相应地computePhysicsResidual中的m,c,k改为从net.Learnables中提取并用exp()解包确保物理量为正m exp(extractLearnable(net.Learnables, m)); c exp(extractLearnable(net.Learnables, c)); k exp(extractLearnable(net.Learnables, k));注意extractLearnable是 MATLAB R2023a 新增函数用于安全获取指定名称的可学习参数。若使用旧版本需遍历net.Learnables结构体查找Parameter字段。3.2 设计多任务损失联合优化位移拟合与参数估计此时损失函数需同时惩罚位移预测误差和参数漂移如先验知识认为 $ k $ 应在 [5,15] 之间。添加 L2 正则项% 提取当前参数 m_est exp(extractLearnable(net.Learnables, m)); c_est exp(extractLearnable(net.Learnables, c)); k_est exp(extractLearnable(net.Learnables, k)); % 参数正则项可选基于先验 reg_m (m_est - 1.0)^2; reg_k (k_est - 10.0)^2; loss_total lambda_d * loss_data lambda_p * loss_physics ... 0.01 * (reg_m reg_k);训练完成后用extractLearnable提取最终参数m_final exp(extractLearnable(net.Learnables, m)); c_final exp(extractLearnable(net.Learnables, c)); k_final exp(extractLearnable(net.Learnables, k)); fprintf(Identified: m%.3f, c%.3f, k%.3f\n, m_final, c_final, k_final);3.3 验证辨识结果用传统 ODE 求解器交叉检验PINN 辨识出的参数必须能复现原始系统行为。用ode45求解经典 ODE并与 PINN 预测对比% 定义 ODE 函数 odeFun (t,x) [x(2); (-c_final*x(2) - k_final*x(1) sin(2*pi*t*0.5))/m_final]; [t_ode, x_ode] ode45(odeFun, t_span, [0; 0]); % 初始位移/速度为0 % PINN 预测需将 t_span 归一化 t_span_norm 2*(t_span - t_min)/(t_max - t_min) - 1; t_span_dl dlarray(t_span_norm, CB); x_pinn double(predict(net, t_span_dl)); % 绘图对比 figure; plot(t_ode, x_ode(:,1), b-, LineWidth,1.5); hold on; plot(t_span, x_pinn, r--, LineWidth,1.5); xlabel(Time (s)); ylabel(Displacement (m)); legend(ODE45, PINN, Location,best); title(sprintf(Parameter ID: m%.2f, c%.2f, k%.2f, m_final, c_final, k_final));若两条曲线高度重合RMSE 1e-3说明 PINN 成功提取了物理本质若存在系统性偏差则需检查物理残差采样密度或增加t_physics点数。4. 提升求解精度与鲁棒性的 4 个关键实践技巧4.1 物理残差采样策略避免频谱泄露与边界奇点物理残差点t_physics的分布直接影响训练稳定性。均匀采样在高频激励下易遗漏关键相位而随机采样可能导致边界区域点稀疏。推荐分段采样法在 $ t \in [0, T/4] $ 和 $ [3T/4, T] $ 区间激励起始/结束用高密度网格步长 0.001在中间区间用低密度网格步长 0.01额外添加 10 个拉丁超立方LHS随机点覆盖整个域。MATLAB 实现t_lhs lhsdesign(1, 10, MaxIterations, 1000); % 10个LHS点 t_lhs t_min (t_max - t_min) * t_lhs; % 映射到实际时间域 t_physics [linspace(0, t_max/4, 50), linspace(3*t_max/4, t_max, 50), ... linspace(t_max/4, 3*t_max/4, 100), t_lhs]; t_physics unique(t_physics); % 去重提示lhsdesign属于 Statistics and Machine Learning Toolbox。若无此工具箱可用rand生成均匀随机点但需确保数量足够≥300。4.2 网络输出后处理强制满足初始条件PINN 网络输出 $ \hat{x}(t) $ 可能不严格满足 $ \hat{x}(0)x_0 $、$ \dot{\hat{x}}(0)v_0 $尤其在数据稀疏时。引入硬约束修正构造辅助函数 $ x_{\text{phys}}(t) x_0 t v_0 t^2 \cdot \hat{x}{\text{net}}(t) $。这样$ x{\text{phys}}(0)x_0 $$ \dot{x}_{\text{phys}}(0)v_0 $ 自动成立且二阶导数仍由网络主导。修改前向传播x_net forward(net, t_dl); x_phys x0 t_dl * v0 t_dl.^2 .* x_net; % t_dl 为归一化时间需映射回实际 t注意t_dl是归一化时间计算t_dl.^2前需还原为实际时间尺度或统一在归一化域内定义初始条件。4.3 损失函数敏感度分析定位主导误差源当训练停滞时需判断是数据拟合不足还是物理约束过强。计算各损失项的梯度范数[~, ~, state] dlfeval(lossFunction, net, t_data_dl, x_data_dl, t_physics_dl, f_physics); grad_data norm(extractGradient(state, output)); % 数据损失梯度 grad_physics norm(extractGradient(state, R_physics)); % 物理损失梯度 fprintf(Gradient norm: Data%.2e, Physics%.2e\n, grad_data, grad_physics);若grad_physics grad_data说明物理残差主导优化应降低lambda_p或增加t_physics点数反之则加强物理约束。4.4 导出为独立函数脱离 Deep Learning Toolbox 运行训练好的 PINN 可导出为纯数值函数供 Simulink 或嵌入式部署使用。利用generateFunctionR2023b或手动提取权重% 提取权重矩阵假设第一层 W1 net.Learnables(1).Value; b1 net.Learnables(2).Value; % 编写纯 MATLAB 函数无 toolbox 依赖 function x_pred pinn_predict(t, W1, b1, W2, b2, W3, b3, W4, b4) t_norm 2*(t - t_min)/(t_max - t_min) - 1; h1 tanh(W1 * t_norm b1); h2 tanh(W2 * h1 b2); h3 tanh(W3 * h2 b3); x_pred W4 * h3 b4; end导出后pinn_predict可在任何 MATLAB 环境包括无 Deep Learning Toolbox 的机器中调用实现零依赖部署。技巧适用场景关键参数/命令效果验证方法分段物理采样高频激励或瞬态响应linspace,lhsdesign物理残差 RMSE 下降 ≥30%初始条件硬约束初始位移/速度已知修改网络输出为 $ x_0 t v_0 t^2 \hat{x}_{\text{net}} $$ x(0) $ 和 $ \dot{x}(0) $ 误差 1e-6梯度范数监控训练发散或停滞extractGradient,norm梯度比值grad_physics/grad_data落入 [0.5, 2.0]权重导出函数部署到无 Toolbox 环境net.Learnables(i).Value, 手写前向传播与predict(net, t_dl)输出差异 1e-10执行完上述步骤你得到的不再是一个黑箱神经网络而是一个严格受牛顿定律约束、参数可解释、结果可验证、部署无依赖的质量-弹簧-阻尼器系统求解器。下一步可将此框架扩展至多自由度系统堆叠网络或图神经网络或耦合热传导、流体阻力等更复杂物理场。本文还有配套的精品资源点击获取