登录论坛

查看完整版本 : [MATLAB工具箱] 求助 fmincon


tjzxd
2009-10-20, 10:58
这是我的语句 总是无法通过,求助哪位高人给指点
function BatchCryOptimzation %间歇结晶
clear all;clc
b=1.45;
kb=285;
Eb=7517;
U=18000;
Hc=44.5;
M=27;
g=1.5;
kg=1.44e+8;
A=0.25;
Eg=4859;
cp=3.8;
pc=2.66e-12;
kv=1.5;
R=8.314;
us0=400;
tf=30;
C0=0.1127;
T0=323;%参数
Tj=100;%降温模型
Tb=[303 323];%温度边界
Cs=.0776+.00246.*Tb-0.0000081.*Tb.^2;
Cm=00.0629+0.00246.*Tb-0.00000714.*Tb.^2;%介稳区
lb=[0 0 0 0 0 303];
ub=[+inf +inf +inf +inf C0 323];
y0=[400 6e+4 2e+7 9e+9 C0 T0];
[y f]=fmincon(@ObjFunc4Fmincon,y0,[],[],[],[],lb,ub,@Nonlcon,[],b,kb,Eb,U,Hc,M,g,kg,A,Eg,cp,pc,kv,R,Tj,tf)%用Fmincon优化
[t,y]=ode45(@Equations,[0 tf],y0,[],b,kb,Eb,U,Hc,M,g,kg,A,Eg,cp,pc,kv,R,Tj);
disp(' t u0 u1 u2 u3 C T '),disp([t,y])
subplot(2,2,1);plot(t,log(y(:,1)),'k--',t,log(y(:,2)),'b:',t,log(y(:,3)),'r-',t,log(y(:,4))),legend('u0','u1','u2','u3'),
xlabel('Time/min'),ylabel('log(ith moment)')
subplot(2,2,2);plot(t,y(:,5)),xlabel('Time/min'),ylabel('Concentration/kgsolute/kg solution')
subplot(2,2,3);plot(t,y(:,6)),xlabel('Time/min'),ylabel('Temperature/K')
subplot(2,2,4);
%.............................
function f=ObjFunc4Fmincon(y0,b,kb,Eb,U,Hc,M,g,kg,A,Eg,cp,pc,kv,R,Tj,tf)
[t,y]=ode45(@Equations,[0 tf],y0,[],b,kb,Eb,U,Hc,M,g,kg,A,Eg,cp,pc,kv,R,Tj); %四=五阶龙格-库塔法
f=y(4);
%........................
function dydt=Equations(t,y,b,kb,Eb,U,Hc,M,g,kg,A,Eg,cp,pc,kv,R,Tj)
u0=y(1);
u1=y(2);
u2=y(3);
u3=y(4);
C=y(5);
T=y(6);
Cs=0.0629+0.00246.*T-0.00000714.*T.^2;
S=(C-Cs)./Cs;
Gt=kg.*exp(-Eg./(R.*T)).*S.^g;
Bt=kb.*exp(-Eb./(R.*T)).*S.^b.*u3;
du0dt=Bt;
du1dt=Gt*u0;
du2dt=2*Gt*u1;
du3dt=3*Gt*u2;
dCdt=-3*pc*kv*Gt*u2;
dTdt=-3*Hc/cp*pc*kv*Gt*u2-U*A/M/cp*(T-Tj);
dydt=[u0;u1;u2;u3;C;T];
%.........................
function [c,ceq]=Nonlcon(y)


c(1)=y(5)-(0.0629+0.00246.*y(6)-0.00000714.*y(6).^2);
c(2)=.0776+.00246.*y(6)-0.0000081.*y(6).^2-y(5);
ceq=[];