matlab含风电场的电力系统的动态经济调度程序，实用性强

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %动态最优潮流 %多时段爬坡约束及基于直流潮流的网络约束 function [fval]=dynamic_optimization_power(wind,temp_gen,temp_load,temp_line) %%%%%%%%%%%%%%%%%%%%%%%输入数据%%%%%%%%%%%%%%%%%%%%%%% num_node=27; slack_id=27; num_sim=1; %模拟次数 delta_t=0.5; temp_wind=wind; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%数据分析%%%%%%%%%%%%%%%%%%%%%%% load_id=temp_load(1,:); temp_load(1,:)=[]; [row_load,column_load]=size(temp_load); num_load=column_load; len=row_load; wind_id=temp_wind(1,:); temp_wind(1,:)=[]; [row_wind,column_wind]=size(temp_wind); num_wind=column_wind; [row_gen,column_gen]=size(temp_gen); num_gen=row_gen; gen_id=temp_gen(:,1); gen_lower=temp_gen(:,2); gen_upper=temp_gen(:,3); gen_ramp=temp_gen(:,4); gen_price=temp_gen(:,5); [row_line,column_line]=size(temp_line); num_line=row_line; line_from=temp_line(:,1)'; line_to=temp_line(:,2)'; line_x=temp_line(:,4)'; line_b=1./line_x; line_limit=temp_line(:,5)'; node_wind=zeros(num_node,num_wind); node_load=zeros(num_node,num_load); for i=1:num_node for j=1:num_wind if wind_id(j)==i node_wind(i,j)=1; end end for k=1:num_load if load_id(k)==i node_load(i,k)=1; end end end node_wind_load=[node_wind,node_load]; node_gen=zeros(num_node,num_gen); for i=1:num_node for j=1:num_gen if gen_id(j)==i node_gen(i,j)=1; end end end line_node=zeros(num_line,num_node); for i=1:num_line from=line_from(i); to=line_to(i); line_node(i,from)=1; line_node(i,to)=-1; end line_node_slackout=line_node; DC_B=zeros(num_node); B_diag=zeros(num_line); for i=1:num_line B_diag(i,i)=1/line_x(i); p=line_from(i); q=line_to(i); DC_B(p,q) = DC_B(p,q)-1/line_x(i); DC_B(q,p) = DC_B(p,q); DC_B(q,q) = DC_B(q,q)+1/line_x(i); DC_B(p,p) = DC_B(p,p)+1/line_x(i); end DC_B0=DC_B; line_node_slackout(:,slack_id)=[]; node_wind_load(slack_id,:)=[]; node_gen(slack_id,:)=[]; DC_B0(slack_id,:) = []; DC_B0(:,slack_id) = []; A0=inv(DC_B0); line_theta=B_diag*line_node_slackout; line_netp=line_theta*A0; line_gen=line_netp*node_gen; %%%%%%%%%%%%%%%%%%%%求最优潮流%%%%%%%%%%%%%%%%%%%% sim_line_cap=zeros(num_line,len); sim_wind_load=zeros(num_node,len); for i=1:len %时段数 for ll=1:num_node for j=1:num_wind if wind_id(j)==ll sim_wind_load(ll,i)=sim_wind_load(ll,i)-temp_wind(i,j); end end for j=1:num_load if load_id(j)==ll sim_wind_load(ll,i)=sim_wind_load(ll,i)+temp_load(i,j); end end end end sim_wind_load(slack_id,:)=[]; for i=1:len %时段数 sim_line_cap(:,i)=line_netp*sim_wind_load(:,i); end Aeq0=zeros(len*num_line+len,len*(num_gen+num_line)); for i=1:len Aeq0((i-1)*num_line+1:i*num_line,(i-1)*(num_gen+num_line)+1:i*(num_gen+num_line))=[line_gen,-eye(num_line)]; end for i=1:len Aeq0(len*num_line+i,(i-1)*(num_gen+num_line)+1:i*(num_gen+num_line))=[ones(1,num_gen),zeros(1,num_line)]; end Beq0=[sim_line_cap(:,1)]; for i=2:len Beq0=[Beq0;sim_line_cap(:,i)]; end for i=1:len Beq0=[Beq0;sum(sim_wind_load(:,i))]; end Lb0=[gen_lower',-line_limit]; Ub0=[gen_upper',line_limit]; for i=2:len Lb0=[Lb0,gen_lower',-line_limit]; Ub0=[Ub0,gen_upper',line_limit]; end Auneq0=zeros(2*(len-1)*num_gen,len*(num_line+num_gen)); for i=1:len-1 Auneq0((i-1)*num_gen+1:i*num_gen,(i-1)*(num_line+num_gen)+1:(i-1)*(num_line+num_gen)+num_gen)=eye(num_gen); Auneq0((i-1)*num_gen+1:i*num_gen,i*(num_line+num_gen)+1:i*(num_line+num_gen)+num_gen)=-eye(num_gen); end for i=1:len-1 Auneq0((len-1)*num_gen+(i-1)*num_gen+1:(len-1)*num_gen+i*num_gen,(i-1)*(num_line+num_gen)+1:(i-1)*(num_line+num_gen)+num_gen)=-eye(num_gen); Auneq0((len-1)*num_gen+(i-1)*num_gen+1:(len-1)*num_gen+i*num_gen,i*(num_line+num_gen)+1:i*(num_line+num_gen)+num_gen)=eye(num_gen); end Buneq0=[gen_ramp'*delta_t]; for i=1:2*len-3 Buneq0=[Buneq0,gen_ramp'*delta_t]; end %此处注释部分为目标函数为二次的求解，而没有注释部分为线性模型，将线性模型注释而将此处二次模型的注释去掉可求解二次问题 % temp=[gen_upper'.^-2,line_limit.^-2]; % for i=1:len-1 % temp=[temp,gen_upper'.^2,line_limit.^2] % end % H=zeros(len*(num_gen+num_line),len*(num_gen+num_line));%*diag(temp) % for i=1:len % H((i-1)*(num_gen+num_line)+1:(i-1)*(num_gen+num_line)+num_gen,(i-1)*(num_gen+num_line)+1:(i-1)*(num_gen+num_line)+num_gen)=2*diag(gen_price_a); % end % C=[gen_price_b',zeros(1,num_line)];%*(gen_upper'.^-2). % for i=1:len-1 % C=[C,gen_price_b',zeros(1,num_line)];%*(gen_upper'.^-2). % end % [sim_gen_line,fval,exitflag,output,lambda]=quadprog(H,C,Auneq0,Buneq0,Aeq0,Beq0,Lb0,Ub0); C=[gen_price',zeros(1,num_line)];%*(gen_upper'.^-2). for i=1:len-1 C=[C,gen_price',zeros(1,num_line)];%*(gen_upper'.^-2). end [sim_gen_line,fval,exitflag,output,lambda]=linprog(C,Auneq0,Buneq0,Aeq0,Beq0,Lb0,Ub0); for i=1:len sim_gen(:,i)=sim_gen_line((i-1)*(num_gen+num_line)+1:(i-1)*(num_gen+num_line)+num_gen); for j=1:num_gen if sim_gen(j,i)<1e-5 sim_gen(j,i)=0; end end end for i=1:len sim_line(:,i)=sim_gen_line((i-1)*(num_gen+num_line)+num_gen+1:i*(num_gen+num_line)); for j=1:num_line if abs(sim_line(j,i))<1e-5 sim_line(j,i)=0; end end end sim_gen=sim_gen'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%