dengchenxy
2010-07-23, 15:16
请假各位高手，我的代码的第二个函数，就是那个二元函数g，目前的求解是没有任何约束的，我要在进行求解时，加入约束t1和t2都大于零，并且t1小于t2来进行求解，之前，我参照了一些教材，但是结果错的很离谱，希望各位高手指点，该如何做。
另外还有一个问题。t0设置成其他参数，就报错了，请问它的设置有什么讲究。
a=3000;
A=50;
c=40;
p=50;
h=3;
b=20;
ie=0.02;
ip=0.04;
B=10;
vb=0:0.05:0.7;
i=1;
diff = 0;
result = zeros(1,15);
ans = zeros(2,15);
%f=@(x)m*(x^2)+n*x+a*c+A/x;
%g=@(t)m1*(t(1)^2)+n1*t(1)+B*b*t(2)+a*B+c*a+c*b*t(2)+m2*(t(1)^2)/t(2)+A/t(2)-a*B*t(1)/t(2)+n2*(t(1)^3)/t(2);
for vb = 0:0.05:0.7;
m = 2/3*b*h+c*b*(1.-vb)*ip-1/3*p*b*ie;
n = b*c+1/2*a*h-1/2*p*a*ie+c*a*(1.-vb)*ip;
m1 = (p*b*ie.*((vb.^2)-1)+p*b*ie*(c/p-(vb.^2)));
n1 = (p*a*ie.*(-1.+vb)+p*a*ie.*(c/p-vb));
m2 = (a*ie.*vb.*(p.*vb-c)-(a*p*ie.*((vb.^2)-(0.5)))+((0.5*a*h)-(b*B)));
n2 = (a*ie.*vb.*(p.*(vb.^2)-c)-p*b*ie.*((vb.^3)-(2/3))+2/3*b*h);
f=@(x)m*(x^2)+n*x+a*c+A/x;
[x(i),y(i)]=fminbnd(f,0,5);
g=@(t)m1*(t(1)^2)+n1*t(1)+B*b*t(2)+a*B+c*a+c*b*t(2)+m2*(t(1)^2)/t(2)+A/t(2)-a*B*t(1)/t(2)+n2*(t(1)^3)/t(2);
t0 = [-1,1];
%options = optimset('GradObj','on');
[t,fval]=fminunc(g,t0);
result(i)=fval;
ans(1,i)=t(1);
ans(2,i)=t(2);
i=i+1;
end
ans,x,y,fval
vb = 0:0.05:0.7;
plot(vb,y)
hold on
vb=0:0.05:0.7;
result;
plot(vb,result)
另外还有一个问题。t0设置成其他参数，就报错了，请问它的设置有什么讲究。
a=3000;
A=50;
c=40;
p=50;
h=3;
b=20;
ie=0.02;
ip=0.04;
B=10;
vb=0:0.05:0.7;
i=1;
diff = 0;
result = zeros(1,15);
ans = zeros(2,15);
%f=@(x)m*(x^2)+n*x+a*c+A/x;
%g=@(t)m1*(t(1)^2)+n1*t(1)+B*b*t(2)+a*B+c*a+c*b*t(2)+m2*(t(1)^2)/t(2)+A/t(2)-a*B*t(1)/t(2)+n2*(t(1)^3)/t(2);
for vb = 0:0.05:0.7;
m = 2/3*b*h+c*b*(1.-vb)*ip-1/3*p*b*ie;
n = b*c+1/2*a*h-1/2*p*a*ie+c*a*(1.-vb)*ip;
m1 = (p*b*ie.*((vb.^2)-1)+p*b*ie*(c/p-(vb.^2)));
n1 = (p*a*ie.*(-1.+vb)+p*a*ie.*(c/p-vb));
m2 = (a*ie.*vb.*(p.*vb-c)-(a*p*ie.*((vb.^2)-(0.5)))+((0.5*a*h)-(b*B)));
n2 = (a*ie.*vb.*(p.*(vb.^2)-c)-p*b*ie.*((vb.^3)-(2/3))+2/3*b*h);
f=@(x)m*(x^2)+n*x+a*c+A/x;
[x(i),y(i)]=fminbnd(f,0,5);
g=@(t)m1*(t(1)^2)+n1*t(1)+B*b*t(2)+a*B+c*a+c*b*t(2)+m2*(t(1)^2)/t(2)+A/t(2)-a*B*t(1)/t(2)+n2*(t(1)^3)/t(2);
t0 = [-1,1];
%options = optimset('GradObj','on');
[t,fval]=fminunc(g,t0);
result(i)=fval;
ans(1,i)=t(1);
ans(2,i)=t(2);
i=i+1;
end
ans,x,y,fval
vb = 0:0.05:0.7;
plot(vb,y)
hold on
vb=0:0.05:0.7;
result;
plot(vb,result)