%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PINN_RS_Main.m
% 一键执行:
% 1. 生成虚拟超声滚压残余应力数据 RS_dataset.csv
% 2. 训练基于物理约束的PINN网络预测残余应力
%
% 作者:GPT-5 (2025)
% MATLAB 2022b+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc; close all;
%% ==================== PART 1 生成虚拟数据 ====================
disp('=== [1/4] Generating synthetic RS dataset ...');
pressureRange = [300 1500]; % N
feedRateRange = [0.06 0.36]; % mm/r
speedRange = [50 300]; % r/min
passRange = [1 6]; % 次
depthPoints = 50; % 每组深度点
numSamples = 90;
rng(1);
dataAll = [];
for i = 1:numSamples
P = rand_range(pressureRange);
F = rand_range(feedRateRange);
S = rand_range(speedRange);
Np = rand_range(passRange);
% 工艺-参数关联关系
A = 200 + 0.1*P - 80*F - 0.1*S + 20*Np;
B = 0.15 + 0.0001*P;
C = 100 + 0.05*P;
D = 3 + 0.002*S;
z0 = 0.25 + 0.05*F;
z = linspace(0,1,depthPoints)';
sigma = -A*exp(-((z-z0).^2)/(2*B^2)) + C*exp(-D*z) + randn(size(z))*20;
for k = 1:length(z)
dataAll = [dataAll; P F S Np sigma(k)];
end
end
header = {'Pressure','FeedRate','Speed','Passes','ResidualStress'};
T = array2table(dataAll,'VariableNames',header);
writetable(T,'RS_dataset.csv');
fprintf('✅ RS_dataset.csv generated (%d samples)\n',height(T));
%% ==================== PART 2 数据预处理 ====================
disp('=== [2/4] Cleaning and preparing dataset ...');
X = T{:,1:4}; Y = T{:,5}; z = linspace(0,1,length(Y))';
% --- LOF异常检测 ---
k = 5;
[~,D] = knnsearch(X,X,'K',k+1);
kDist = D(:,end);
reachDist = max(kDist(D(:,end)), D);
LRD = 1 ./ mean(reachDist,2);
LOF = mean(LRD(knnsearch(X,X,'K',k+1)),2) ./ LRD;
mask = LOF < 1;
X = X(mask,:); Y = Y(mask);
Xn = normalize(X);
Yn = normalize(Y);
fprintf('✅ %d valid samples after cleaning\n',length(Y));
%% ==================== PART 3 定义PINN模型 ====================
disp('=== [3/4] Building and training PINN model ...');
% 定义组合函数与导数
sigma_combined = @(z,A,B,C,D,z0) ...
-A.*exp(-((z-z0).^2)./(2*B.^2)) + C.*exp(-D.*z);
dSigma = @(z,A,B,C,D,z0) ...
A.*(z-z0)./B.^2 .* exp(-((z-z0).^2)/(2*B.^2)) - C.*D.*exp(-D.*z);
bfgs_maxsigma = @(A,B,C,D,z0) bfgs_find_max(A,B,C,D,z0,dSigma);
% 网络结构 (简化Structure 1+2)
inputSize = 4; hiddenSize = 50;
layers = [
featureInputLayer(inputSize,"Name","input")
fullyConnectedLayer(hiddenSize)
reluLayer
fullyConnectedLayer(hiddenSize)
reluLayer
fullyConnectedLayer(1,"Name","output")];
lgraph = layerGraph(layers);
net = dlnetwork(lgraph);
% 训练参数
learnRate = 5e-4; numEpochs = 300; miniBatchSize = 32;
gradDecay = 0.9; sqGradDecay = 0.999;
trailingAvg = []; trailingAvgSq = [];
numIterPerEpoch = floor(size(Xn,1)/miniBatchSize);
for epoch = 1:numEpochs
idx = randperm(size(Xn,1));
for i = 1:numIterPerEpoch
ind = idx((i-1)*miniBatchSize+1 : i*miniBatchSize);
XBatch = dlarray(Xn(ind,:)','CB');
YBatch = dlarray(Yn(ind,:)','CB');
[loss,gradients] = dlfeval(@modelLoss,net,XBatch,YBatch,z,...
sigma_combined,dSigma,bfgs_maxsigma);
[net,trailingAvg,trailingAvgSq] = adamupdate(net,gradients,...
trailingAvg,trailingAvgSq,epoch,learnRate,gradDecay,sqGradDecay);
end
if mod(epoch,10)==0
fprintf("Epoch %d/%d, Loss = %.5f\n",epoch,numEpochs,double(extractdata(loss)));
end
end
%% ==================== PART 4 模型评估与可视化 ====================
disp('=== [4/4] Evaluating and visualizing results ...');
Y_pred = predict(net,Xn')';
Y_pred = rescale(Y_pred, min(Y), max(Y));
RMSE = sqrt(mean((Y - Y_pred).^2));
R2 = corr(Y,Y_pred)^2;
fprintf("\n✅ PINN Performance: R² = %.3f, RMSE = %.3f MPa\n",R2,RMSE);
figure;
scatter(Y,Y_pred,15,'filled'); grid on;
xlabel('Measured Residual Stress (MPa)');
ylabel('Predicted Residual Stress (MPa)');
title(sprintf('PINN Residual Stress Prediction\nR²=%.3f, RMSE=%.1f MPa',R2,RMSE));
figure;
subplot(1,2,1);
plot(Y(1:500),'k','LineWidth',1.2); hold on;
plot(Y_pred(1:500),'r--','LineWidth',1.2);
xlabel('Sample Index'); ylabel('Residual Stress (MPa)');
legend('Measured','Predicted'); grid on;
title('Comparison (Partial Samples)');
subplot(1,2,2);
err = Y - Y_pred;
histogram(err,30);
xlabel('Prediction Error (MPa)');
ylabel('Count'); grid on;
title('Error Distribution');
disp('✅ All steps completed successfully.');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ===================== 函数定义部分 =====================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [loss,gradients] = modelLoss(net,XBatch,YBatch,z,sigma_combined,dSigma,bfgs_maxsigma)
preds = forward(net,XBatch);
preds = preds(1,:);
A=1;B=0.2;C=0.5;D=0.3;z0=0.2; % 简化参数
zmax = bfgs_maxsigma(A,B,C,D,z0);
delta = 1;
err = preds - YBatch;
lossData = mean((abs(err)<delta).*0.5.*err.^2 + (abs(err)>=delta).*delta.*(abs(err)-0.5*delta));
dSig = dSigma(z,A,B,C,D,z0);
F1 = max(0,dSig(z<zmax));
F2 = max(0,-dSig(z>zmax));
lossPhys = mean([F1;F2].^2);
loss = lossData + 0.1*lossPhys;
gradients = dlgradient(loss,net.Learnables);
end
function zmax = bfgs_find_max(A,B,C,D,z0,dSigma)
z = z0; Bk=1; tol=1e-6; maxIter=100;
for k=1:maxIter
grad = dSigma(z,A,B,C,D,z0);
if abs(grad)<tol, break; end
p = -Bk*grad;
alpha = 0.1;
z_new = z + alpha*p;
yk = dSigma(z_new,A,B,C,D,z0)-grad;
sk = z_new - z;
Bk = Bk + (yk^2)/(yk*sk) - (Bk^2*sk^2)/(Bk*sk^2);
z = z_new;
end
zmax = z;
end
function val = rand_range(range)
val = range(1) + (range(2)-range(1))*rand;
end
修改上述代码,使matlab能够运行