[求助]那位大侠能知道小弟程序那里有问题？

[jackwen]

那位大侠能帮助小弟，该程序，利用fmincon解决电网中发电节点的问题，我从30节点改到9节点，还是出问题，但问题是同一个问题，快要崩溃了，不知道怎么解决。
问题：
Warning: Large-scale (trust region) method does not currently solve this type of problem,
switching to medium-scale (line search).
> In fmincon at 260
In robus at 45
这是警告。
x =

Columns 1 through 12

50.6789 31.4128 22.7126 -0.6789 2.4679 0.2478 0.6789 -2.4679 -0.2478 9.4371 -0.5761 -0.9122

Columns 13 through 18

1.2854 -2.1539 1.4942 -2.1043 1.4744 -2.3799


fval =

428.9764


exitflag =

0


p =

iterations: 31
funcCount: 1019
stepsize: 0.0625
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
firstorderopt: 13.2031
cgiterations: []
message: [1x79 char]
这是结果，不收敛。exitflag =0

下面是我的原程序：
第一个，目标函数function y=myfun(x)
%% generator cost data
% 1 startup shutdown n x0 y0 ... xn yn
% 2 startup shutdown n c(n-1) ... c0
gencost =[
2 1500 0 3 0.11 5 150;
2 2000 0 3 0.085 1.2 600;
2 3000 0 3 0.1225 1 335;
];
n1=0.001;
n2=-n1;
M=gencost(:,4);
a=gencost(:,3);
b=gencost(:,6);
c=gencost(:,5);
y1=a(1)+b(1)*x(1)+c(1)*x(1)^2+a(2)+b(2)*x(2)+c(2)*x(2)^2+a(3)+b(3)*x(3)+c(3)*x(3)^2;
y2=M(1)*x(1)+M(2)*x(2)+M(3)*x(3);
y3=x(4)*n1*M(1)+x(5)*n1*M(2)+x(6)*n1*M(3);
y4=x(7)*n2*M(1)+x(8)*n2*M(2)+x(9)*n2*M(3);
y=y1-y2-y3+y4;

第二个约束程序
function[c,ceq]=mycon(x)
%% system MVA base
baseMVA = 100;
%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
bus = [
1 3 0 0 0 0 1 1 0 345 1 1.1 0.9;
2 2 0 0 0 0 1 1 0 345 1 1.1 0.9;
3 2 0 0 0 0 1 1 0 345 1 1.1 0.9;
4 1 0 0 0 0 1 1 0 345 1 1.1 0.9;
5 1 90 30 0 0 1 1 0 345 1 1.1 0.9;
6 1 0 0 0 0 1 1 0 345 1 1.1 0.9;
7 1 100 35 0 0 1 1 0 345 1 1.1 0.9;
8 1 0 0 0 0 1 1 0 345 1 1.1 0.9;
9 1 125 50 0 0 1 1 0 345 1 1.1 0.9;
];
%% generator data
% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin
gen = [
1 0 0 300 -300 1 100 1 250 10;
2 163 0 300 -300 1 100 1 300 10;
3 85 0 300 -300 1 100 1 270 10;
];

%% branch data
% fbus tbus r x b rateA rateB rateC ratio angle status
branch = [
1 4 0 0.0576 0 250 250 250 0 0 1;
4 5 0.017 0.092 0.158 250 250 250 0 0 1;
5 6 0.039 0.17 0.358 150 150 150 0 0 1;
3 6 0 0.0586 0 300 300 300 0 0 1;
6 7 0.0119 0.1008 0.209 150 150 150 0 0 1;
7 8 0.0085 0.072 0.149 250 250 250 0 0 1;
8 2 0 0.0625 0 250 250 250 0 0 1;
8 9 0.032 0.161 0.306 250 250 250 0 0 1;
9 4 0.01 0.085 0.176 250 250 250 0 0 1;
];
pd=bus(:,3);
pg=[x(1) x(2) x(3) 0 0 0 0 0 0];%9个数据
g=0;
k=0;
l=0;
lb2=[50 20 15];
ub2=[200 80 50];
T_max=branch(:,6);
for i=1:3
ceq(i)=x(i)-x(10)-x(i+3)+x(i+6);
c(i)=-x(i+3);
c(i+3)=-x(i+6);
end
%调用函数makeBdc计算直流导纳矩阵
B=makeBdc(baseMVA,bus,branch);
%直流潮流约束方程
Num=9;
for i=2:Num
for j=2:Num
if j==i
h=B(i,j)*x(i+9);
else
h=(-1)*B(i,j)*x(j+9);
end
g=g+h;
h=0;
end
ceq(i+2)=g+pd(i)-pg(i);
g=0;
end
%编写DC-network model约束，具体从c(7)开始 Num=9(修改此约束12.4)
for i=1:Num
if i==1
Fa=0;
else
Fa=1;
end
for j=1:Num
if j==1
Fb=0;
else
Fb=1;
end
h=B(i,j)*(Fa*x(i+9)-Fb*x(j+9));
k=k+h;
h=0;
end
c(i+6)=-(k+pg(i));
c(i+15)=k+pg(i)-pd(i);
k=0;
end
%编写Intact-network line-flow constraints约束，具体从c()开始 40个a(i,j)
for i=1:Num
if i==1
Fa=0;
else
Fa=1;
end
for j=1:Num
if j==1
Fb=0;
else
Fb=1;
end
c(i+24)=B(i,j)*(Fa*x(i+9)-Fb*x(j+9))-T_max(i);
c(i+33)=-B(i,j)*(Fa*x(i+9)-Fb*x(j+9))-T_max(i);
end
end
for i=1:3
c(i+42)=x(i)-ub2(i);
c(i+45)=lb2(i)-x(i);
end
return
第三个 执行程序
function [x,fval,exitflag,p]=robus
lb1=zeros(18,1);
ub1=ones(18,1);
lb2=[50 20 15];
ub2=[200 80 50];
for i=4:18
lb2(i)=0;
ub2(i)=0;
end
%发电成本参数
gencost =[
2 1500 0 3 0.11 5 150;
2 2000 0 3 0.085 1.2 600;
2 3000 0 3 0.1225 1 335;
];
n1=0.001;
n2=-n1;
for i=1:3
e(i)=lb2(i);
end
for i=4:10
e(i)=0;
end
%% bus data
% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin
bus = [
1 3 0 0 0 0 1 1 0 345 1 1.1 0.9;
2 2 0 0 0 0 1 1 0 345 1 1.1 0.9;
3 2 0 0 0 0 1 1 0 345 1 1.1 0.9;
4 1 0 0 0 0 1 1 0 345 1 1.1 0.9;
5 1 90 30 0 0 1 1 0 345 1 1.1 0.9;
6 1 0 0 0 0 1 1 0 345 1 1.1 0.9;
7 1 100 35 0 0 1 1 0 345 1 1.1 0.9;
8 1 0 0 0 0 1 1 0 345 1 1.1 0.9;
9 1 125 50 0 0 1 1 0 345 1 1.1 0.9;
];
for i=11:18
e(i)=bus(i-10,9);
%e(i)=0;
end
x0=e;
p=e;
%调用fmincon函数
options=optimset('Display','iter','FunValCheck','on','MaxFunEvals',1000,'MaxIter',1000,'TolFun',1e-1,'TolX',1e-1);
[x,fval,exitflag,p]=fmincon(@myfun,x0,[],[],[],[],[],[],@mycon,options);
那位大侠能帮助小弟看看，不胜感激。