应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 缩聚反应动力学参数估计 置信区间太大 求程序改进
51/1返回列表
查看: 774  |  回复: 2
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

cyzeng

金虫 (小有名气)

[求助] 缩聚反应动力学参数估计 置信区间太大 求程序改进

缩聚反应动力学参数估计 置信区间太大 求程序改进
程序
clc;clear;
format long
global K1 K2 K3 cin rin T R
T=[483.15 493.15 503.15 513.15 523.15];%动力学实验温度点，K
R=8.3145;%J/(K*mol)
lb=[1e2   3e4  1e2   3e4 1e2   3e4];
ub=[1e10  8e4  1e10  8e4 1e10  8e4];

%实验测定平衡常数随温度变化
deltS1=43.2237;
deltH1=1.58e+4;
deltS2=50.4599;
deltH2=2.01e+4;
K1=exp(deltS1/R-deltH1./(R*T));
K2=exp(deltS2/R-deltH2./(R*T));
K3=K2./K1;
% 动力学数据-不同温度下间歇反应器内各浓度随时间变化
%           time/h    c1       c2      c3      c4      c5
ExpData_210=[0      0.372         0.076         1.630         4.439         0.0249
             1      0.146         0.087         1.339         4.718         0.0127
             2      0.059         0.073         1.270         4.820         0.0100
             3      0.025         0.070         1.245         4.856         0.0052];
ExpData_220=[0      0.400         0.079         1.683         4.387         0.0306;
             1      0.125         0.054         1.016         4.990         0.0106;
             2      0.042         0.033         0.948         5.100         0.0067;
             3      0.015         0.028         0.886         5.162        0.0033];
ExpData_230=[0      0.290         0.070         1.620         4.488         0.0239
             0.5    0.102         0.061         1.071         4.954         0.0099
             1      0.040         0.039         1.033         5.035         0.0054
             2      0.008         0.036         1.011         5.385         0.0030];
ExpData_240=[0      0.285         0.065         1.386         4.656         0.0229
             0.5    0.088         0.033         0.823         5.165         0.0086
             1      0.035         0.028         0.781         5.224         0.0043
             1.5    0.016         0.025         0.750         5.258         0.0020];
ExpData_250=[0      0.228         0.050         1.090         4.902         0.0236
             0.5    0.054         0.018         0.635         5.329         0.0086
             0.75   0.038   0.016   0.608   5.357   0.0064
             1      0.024         0.014         0.589         5.378         0.0041];           
%生成动力学数据三维矩阵-转换成行向量
ExpData=ones(6,4,5);
ExpData(:,:,1)=ExpData_210';
ExpData(:,:,2)=ExpData_220';
ExpData(:,:,3)=ExpData_230';
ExpData(:,:,4)=ExpData_240';
ExpData(:,:,5)=ExpData_250';
time=ExpData(1,:,;
cexp=ExpData(2:6,:,;
c0=cexp(:,1,;        
%浓度插值及微分求导得反应速率
for i=1:5
    [cin(:,:,i) rin(:,:,i)]=Rate(time(:,:,i),cexp(:,:,i));
end
cin;
rin;
%拟合EG和H2O反应速率曲线
for i=1:5
    ti2=linspace(time(:,1,i),time(:,4,i),61);
    ppp2(:,:,i)=polyfit(ti2,rin(2,:,i),4);
    ppp5(:,:,i)=polyfit(ti2,rin(5,:,i),4);
end
ppp2;ppp5;

b0=[2000.4 30547.1 16289.2 42486.9 13232.6 48038.2];

%使用fmincon进行参数估计
% options=optimset('MaxFunEvals',2400);
[b,fval,flag]=fmincon(@objfuncfmincon,b0,[],[],[],[],lb,ub,[],[],c0,cin,ppp2,ppp5);
fval
b_fmincon=b;
b0=b_fmincon
%以fmincon估计的k0,使用lsqnonlin进行参数估计
options=optimset('MaxFunEvals',2400);
[b,resnorm,resid,exitflag,output,lambda,jacobian]=lsqnonlin(@objfunclsqnonlin,b0,lb,ub,options,c0,cin,ppp2,ppp5);
resnorm
b
bi=nlparci(b,resid,jacobian);
fprintf('\tk10=%.1f±%.1f\n',b(1),bi(1,2)-b(1))
fprintf('\tEa1=%.1f±%.1f\n',b(2),bi(2,2)-b(2))
fprintf('\tk30=%.1f±%.1f\n',b(3),bi(3,2)-b(3))
fprintf('\tEa3=%.1f±%.1f\n',b(4),bi(4,2)-b(4))
fprintf('\tk50=%.1f±%.1f\n',b(5),bi(5,2)-b(5))
fprintf('\tEa5=%.1f±%.1f\n',b(6),bi(6,2)-b(6))
%模型计算值与实验点比较
%210度
ti3=linspace(time(:,1,1),time(:,4,1),61);
[t,c]=ode45(@masslsq1,ti3,c0(:,:,1),[],b,ppp2,ppp5);
figure
plot(time(:,:,1),cexp(1,:,1),'ks',time(:,:,1),cexp(2,:,1),'k<',time(:,:,1),cexp(3,:,1),'k>',time(:,:,1),cexp(4,:,1),'k*',time(:,:,1),cexp(5,:,1),'ko',...
    ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=483.15K')
%220度
ti3=linspace(time(:,1,2),time(:,4,2),61);
[t,c]=ode45(@masslsq2,ti3,c0(:,:,2),[],b,ppp2,ppp5);
figure
plot(time(:,:,2),cexp(1,:,2),'ks',time(:,:,2),cexp(2,:,2),'k<',time(:,:,2),cexp(3,:,2),'k>',time(:,:,2),cexp(4,:,2),'k*',time(:,:,2),cexp(5,:,2),'ko',...
    ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=493.15K')
%230度
ti3=linspace(time(:,1,3),time(:,4,3),61);
[t,c]=ode45(@masslsq3,ti3,c0(:,:,3),[],b,ppp2,ppp5);
figure
plot(time(:,:,3),cexp(1,:,3),'ks',time(:,:,3),cexp(2,:,3),'k<',time(:,:,3),cexp(3,:,3),'k>',time(:,:,3),cexp(4,:,3),'k*',time(:,:,3),cexp(5,:,3),'ko',...
    ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=503.15K')
%240度
ti3=linspace(time(:,1,4),time(:,4,4),61);
[t,c]=ode45(@masslsq4,ti3,c0(:,:,4),[],b,ppp2,ppp5);
figure
plot(time(:,:,4),cexp(1,:,4),'ks',time(:,:,4),cexp(2,:,4),'k<',time(:,:,4),cexp(3,:,4),'k>',time(:,:,4),cexp(4,:,4),'k*',time(:,:,4),cexp(5,:,4),'ko',...
    ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=513.15K')
%250度
ti3=linspace(time(:,1,5),time(:,4,5),61);
[t,c]=ode45(@masslsq5,ti3,c0(:,:,5),[],b,ppp2,ppp5);
figure
plot(time(:,:,5),cexp(1,:,5),'ks',time(:,:,5),cexp(2,:,5),'k<',time(:,:,5),cexp(3,:,5),'k>',time(:,:,5),cexp(4,:,5),'k*',time(:,:,5),cexp(5,:,5),'ko',...
    ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=523.15K')

%定义插值和反应速率函数
function [cin rin]=Rate(t,c)
% knots=1; K=4;
% sp=spap2(knots,K,t,c);
% pp=fnder(sp);
% ti=linspace(t(1),t(end),61);
% cin=fnval(sp,ti);
% rin=fnval(pp,ti);           
sp=csaps(t,c,0.99);
pp=fnder(sp);
ti=linspace(t(1),t(end),61);
cin=fnval(sp,ti);
rin=fnval(pp,ti);

%定义objfuncfmincon
function ffmi=objfuncfmincon(b,c0,cin,ppp2,ppp5)
tspan=linspace(0,3,61);
[t,c]=ode45(@masslsq1,tspan,c0(:,:,1),[],b,ppp2,ppp5);
f1=sum((c(:,1)-cin(1,:,1)').^2)+sum((c(:,3)-cin(3,:,1)').^2)+sum((c(:,4)-cin(4,:,1)').^2);
[t,c]=ode45(@masslsq2,tspan,c0(:,:,2),[],b,ppp2,ppp5);
f2=sum((c(:,1)-cin(1,:,2)').^2)+sum((c(:,3)-cin(3,:,2)').^2)+sum((c(:,4)-cin(4,:,2)').^2);
tspan2=linspace(0,2,61);
[t,c]=ode45(@masslsq3,tspan2,c0(:,:,3),[],b,ppp2,ppp5);
f3=sum((c(:,1)-cin(1,:,3)').^2)+sum((c(:,3)-cin(3,:,3)').^2)+sum((c(:,4)-cin(4,:,3)').^2);
tspan3=linspace(0,1.5,61);
[t,c]=ode45(@masslsq4,tspan3,c0(:,:,4),[],b,ppp2,ppp5);
f4=sum((c(:,1)-cin(1,:,4)').^2)+sum((c(:,3)-cin(3,:,4)').^2)+sum((c(:,4)-cin(4,:,4)').^2);
tspan4=linspace(0,1,61);
[t,c]=ode45(@masslsq5,tspan4,c0(:,:,5),[],b,ppp2,ppp5);
f5=sum((c(:,1)-cin(1,:,5)').^2)+sum((c(:,3)-cin(3,:,5)').^2)+sum((c(:,4)-cin(4,:,5)').^2);
ffmi=f1+f2+f3+f4+f5;

%定义objfunclsqnonlin
function flsq=objfunclsqnonlin(b,c0,cin,ppp2,ppp5)
tspan=linspace(0,3,61);
[t,c]=ode45(@masslsq1,tspan,c0(:,:,1),[],b,ppp2,ppp5);
f1=c(:,1)-cin(1,:,1)';
f2=c(:,3)-cin(3,:,1)';
f3=c(:,4)-cin(4,:,1)';
[t,c]=ode45(@masslsq2,tspan,c0(:,:,2),[],b,ppp2,ppp5);
f4=c(:,1)-cin(1,:,2)';
f5=c(:,3)-cin(3,:,2)';
f6=c(:,4)-cin(4,:,2)';
tspan2=linspace(0,2,61);
[t,c]=ode45(@masslsq3,tspan2,c0(:,:,3),[],b,ppp2,ppp5);
f7=c(:,1)-cin(1,:,3)';
f8=c(:,3)-cin(3,:,3)';
f9=c(:,4)-cin(4,:,3)';
tspan3=linspace(0,1.5,61);
[t,c]=ode45(@masslsq4,tspan3,c0(:,:,4),[],b,ppp2,ppp5);
f10=c(:,1)-cin(1,:,4)';
f11=c(:,3)-cin(3,:,4)';
f12=c(:,4)-cin(4,:,4)';
tspan4=linspace(0,1,61);
[t,c]=ode45(@masslsq5,tspan4,c0(:,:,5),[],b,ppp2,ppp5);
f13=c(:,1)-cin(1,:,5)';
f14=c(:,3)-cin(3,:,5)';
f15=c(:,4)-cin(4,:,5)';
flsq=[f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15];

%定义masslsq1-210度
function dcdt1=masslsq1(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(1)));
k2=k1/K1(1);
k3=b(3)*exp(-b(4)/(R*T(1)));
k4=k3/K2(1);
k5=b(5)*exp(-b(6)/(R*T(1)));
k6=k5/K3(1);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,1)+ppp2(:,4,1)*t+ppp2(:,3,1)*t.^2+ppp2(:,2,1)*t.^3+ppp2(:,1,1)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,1)+ppp5(:,4,1)*t+ppp5(:,3,1)*t.^2+ppp5(:,2,1)*t.^3+ppp5(:,1,1)*t.^4;
dcdt1=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq2-220度
function dcdt2=masslsq2(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(2)));
k2=k1/K1(2);
k3=b(3)*exp(-b(4)/(R*T(2)));
k4=k3/K2(2);
k5=b(5)*exp(-b(6)/(R*T(2)));
k6=k5/K3(2);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,2)+ppp2(:,4,2)*t+ppp2(:,3,2)*t.^2+ppp2(:,2,2)*t.^3+ppp2(:,1,2)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,2)+ppp5(:,4,2)*t+ppp5(:,3,2)*t.^2+ppp5(:,2,2)*t.^3+ppp5(:,1,2)*t.^4;
dcdt2=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq3-230度
function dcdt3=masslsq3(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(3)));
k2=k1/K1(3);
k3=b(3)*exp(-b(4)/(R*T(3)));
k4=k3/K2(3);
k5=b(5)*exp(-b(6)/(R*T(3)));
k6=k5/K3(3);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,3)+ppp2(:,4,3)*t+ppp2(:,3,3)*t.^2+ppp2(:,2,3)*t.^3+ppp2(:,1,3)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,3)+ppp5(:,4,3)*t+ppp5(:,3,3)*t.^2+ppp5(:,2,3)*t.^3+ppp5(:,1,3)*t.^4;
dcdt3=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq4-240度
function dcdt4=masslsq4(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(4)));
k2=k1/K1(4);
k3=b(3)*exp(-b(4)/(R*T(4)));
k4=k3/K2(4);
k5=b(5)*exp(-b(6)/(R*T(4)));
k6=k5/K3(4);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,4)+ppp2(:,4,4)*t+ppp2(:,3,4)*t.^2+ppp2(:,2,4)*t.^3+ppp2(:,1,4)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,4)+ppp5(:,4,4)*t+ppp5(:,3,4)*t.^2+ppp5(:,2,4)*t.^3+ppp5(:,1,4)*t.^4;
dcdt4=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq5-250度
function dcdt5=masslsq5(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(5)));
k2=k1/K1(5);
k3=b(3)*exp(-b(4)/(R*T(5)));
k4=k3/K2(5);
k5=b(5)*exp(-b(6)/(R*T(5)));
k6=k5/K3(5);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,5)+ppp2(:,4,5)*t+ppp2(:,3,5)*t.^2+ppp2(:,2,5)*t.^3+ppp2(:,1,5)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,5)+ppp5(:,4,5)*t+ppp5(:,3,5)*t.^2+ppp5(:,2,5)*t.^3+ppp5(:,1,5)*t.^4;
dcdt5=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

结果
Optimization terminated: directional derivative predicts change in
  objective value less than options.TolFun and maximum constraint violation
  is less than options.TolCon.
No active inequalities.

fval =

   9.975490368849357


b0 =

  1.0e+004 *

   0.021390412595752   3.140646199334575   1.805623620390919   3.561747395742712   1.562837065436626   4.047673875392536

Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.

resnorm =

   8.289356826509787


b =

  1.0e+007 *

   0.000824325058983   0.003320473746231   0.141719950795143   0.005636884012449   3.314837171169622   0.006838741499027

        k10=8243.3±936367.3
        Ea1=33204.7±461448.7
        k30=1417199.5±20061412.4
        Ea3=56368.8±58203.6
        k50=33148371.7±584453034.8
        Ea5=68387.4±72746.8
回复此楼

» 猜你喜欢
1楼 2012-09-07 10:17:36
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

梅子炒米粉

禁虫
本帖内容被屏蔽
3楼2021-04-24 21:40:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 3 个回答
查看全部 3 个回答
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 085600 英一数二272求调剂 5+6 vida_a 2026-03-01 15/750 2026-03-02 18:50 by caszguilin
[考研] 材料调剂 +3 恒顺自然 2026-03-02 3/150 2026-03-02 18:49 by L135790
[考研] 高分子化学与物理调剂 +6 好好好1233 2026-02-28 15/750 2026-03-02 18:47 by caszguilin
[考研] 283求调剂 +4 鹿沫笙 2026-03-02 4/200 2026-03-02 18:41 by caszguilin
[考研] 考研复试调剂，过国家线的同学都可报名 +4 黑！在干嘛 2026-02-28 5/250 2026-03-02 16:48 by chock1337
[考研] 291 求调剂 +3 化工2026届毕业� 2026-03-02 3/150 2026-03-02 12:55 by houyaoxu
[考研] 295求调剂 +8 19171856320 2026-02-28 8/400 2026-03-02 11:19 by yuchj
[基金申请] 此成果不能导入原因：元数据必填信息不完整，可 进行补充。 +4 Kittylucky 2026-03-02 5/250 2026-03-02 11:07 by jurkat.1640
[考研] 281求调剂 +5 2026计算机_诚心 2026-03-01 8/400 2026-03-02 11:05 by 汪！？！
[考研] 调剂 +3 13853210211 2026-03-02 4/200 2026-03-02 10:16 by 13853210211
[基金申请] 本子写完了，给DS兄弟看了，得了92分 +3 Doma 2026-03-01 7/350 2026-03-02 00:00 by jnzsy
[基金申请] 成果系统访问量大，请一小时后再尝试。---NSFC啥时候好哦，已经两天这样了 +4 NSFC2026我来了 2026-02-28 4/200 2026-03-01 22:37 by 铁门栓
[硕博家园] 博士自荐 +7 科研狗111 2026-02-26 11/550 2026-03-01 22:24 by 哲平L
[考研] 299求调剂 +3 Y墨明棋妙Y 2026-02-28 5/250 2026-03-01 21:01 by tangxiaotian
[考研] 298求调剂 +6 axyz3 2026-02-28 6/300 2026-03-01 19:00 by 18137688336
[考研] 290求调剂 +9 材料专硕调剂； 2026-02-28 11/550 2026-03-01 17:21 by sunny81
[考研] 321求调剂一志愿东北林业大学材料与化工英二数二 +4 虫虫虫虫虫7 2026-03-01 7/350 2026-03-01 16:52 by caszguilin
[考研] 285求调剂 +8 满头大汗的学生 2026-02-28 8/400 2026-03-01 16:47 by caszguilin
[论文投稿] Optics letters投稿被拒求助 30+3 luckyry 2026-02-26 4/200 2026-03-01 09:06 by babero
[硕博家园] 【博士招生】太原理工大学2026化工博士 +4 N1ce_try 2026-02-24 8/400 2026-02-26 08:40 by N1ce_try
信息提示
关闭
请填处理意见
关闭