1.算法描述

高载能企业执行子站接收负荷调整指令后，需将有功功率、无功功率调整总量合理分配给各用能设备/系统。研究高载能负荷响应优化控制模型，建立以高载能企业响应效益最优为目标，以各用能设备/系统在不同时间尺度下有功、无功响应容量为变量，以各用能设备/系统在不同时间尺度下响应潜力、无功补偿容量、用电重要等级、生产工艺环节协调配合、有功功率及无功功率调整总量等为约束，提出寻优求解方法。

对于铝电解厂来说，风电引入到铝电解厂，但是铝电解厂有自备电厂，正常情况下是自备电厂恒定地供给铝电解厂，铝电解厂的功率主要可分为两种设备，一个是主要设备铝电解槽（占总功率的95%），其余是辅助设备消耗的功率（约占总功率的5%）。正常工作条件下，铝电解槽和辅助系统24小时的功率如下面两图所示，基本上处于稳定状态，此时P铝电解槽+P辅助= P自备电厂。

现在如果把风电引入到铝电解厂了，则原来的自备电厂供电的功率会发生变化，同时铝电解厂消耗的总功率也会发生变化。现在:

P铝电解槽+P辅助=P风电+P自备电厂

如果现在预先给定了引入铝电解厂风电24小时的曲线，如下图所示，现在在这24小时期间，如何调整自备电厂的发电、以及如何将自备电厂和风电的总功率分配给主要设备铝电解槽和辅助系统。才能使在引入风电的情况下电解厂获得的效益（赚的钱）最高，这时，最好使电解槽的功率得以调高，多增加铝的产量，会多带来效益。假设风电的价格是0.2元/kWh，自备电厂的价格是0.3元/kWh，铝的价格是13000元/吨。

一些约束条件可以设置如下：

主要设备电解槽电解槽的额定功率是700MW，调整时不能低于730MW，不能高于770MW；

辅助设备的额定功率是 50MW，调整时不能低于45 MW，不能高于55MW；

自备电厂的额定输出功率为750MW，每次调整不能低于600MW，不能高于780MW。

根据上面这幅图给出的风电出力曲线，经过对铝电解厂和发电厂负荷的调整，使得这24小时中盈利最大。并求出发电厂、每种负荷（铝电解槽和辅助设备）24小时的调整曲线。

2.仿真效果预览

matlab2022a仿真如下：

3.MATLAB核心程序


%输入风电变换情况
Times = [0.25:0.25:24];
Pfd0  = [3*ones(1,20),5*ones(1,12),10*ones(1,20),12*ones(1,8),8*ones(1,24),4*ones(1,12)]*1e6;
figure;
plot(Times,Pfd0/1e6,'r','linewidth',2);
xlabel('Time(h)');
ylabel('MW');
axis([0,24,0,15]);

Timesdelay = 2*4;%2小时后参与调解

Pdc21  = [];
Pfl21  = [];
Pdzl21 = [];
Pdc22  = [];
Pfl22  = [];
Pdzl22 = [];


for t = 1:length(Times)

%定义初始值，初始值
Pdc   = 40e6;%电厂调节功率
Pfl   = 15e6;%辅助设备调节功率
Pdzl  = 25e6;%电炉调节功率
Pfd   = Pfd0(t);

%风电价格波动
Pricefd = BD(t); 

%%
%下面开始使用遗传优化算法
%根据遗传算法进行参数的拟合
MAXGEN = 80;
NIND   = 200;
Nums   = 2;

Chrom  = crtbp(NIND,Nums*10);
%各个变量的约束条件,这几个参数稍微改了下，否则加入众多约束后，收敛非常慢
A1     = 20e6;
B1     = 60e6;

A2     = 12e6;
B2     = 20e6;

A3     = 30e6;
B3     = 35e6;
Areas  = [A2, A3;
          B2, B3];

FieldD = [rep([10],[1,Nums]);Areas;rep([0;0;0;0],[1,Nums])];

gen    = 0;
for a=1:1:NIND 
    %计算对应的目标值
    Moneys  = func_obj(Pdc,Pfl,Pdzl,Pfd,t,Pricefd);
    Js(a,1) = Moneys;
end

Objv  = (Js+eps);
gen   = 0; 
Pdc2  = [];
Pfl2  = [];
Pdzl2 = [];
%%
while gen < MAXGEN;   
      gen
      t
      
      flag=0;
      while flag == 0
      Pe0 = 0.8;
      pe1 = 0.001; 
      FitnV=ranking(Objv);    
      Selch=select('sus',Chrom,FitnV);    
      Selch=recombin('xovsp', Selch,Pe0);   
      Selch=mut( Selch,pe1);   
      phen1=bs2rv(Selch,FieldD);   
      NS=0;
          for aa=1:NIND
              if phen1(a,1) +  phen1(a,2) - Pfd <= B1 & phen1(a,1) +  phen1(a,2) - Pfd >= A1
                 NS=NS+1;
              end
          end
          if NS > NIND/2
             flag = 1;
          else
             flag = 0; 
          end
      end
      Pdc   = [];
      Pfl   = [];
      Pdzl  = [];
      for a=1:1:NIND  
          Pfl(a)   = phen1(a,1);
          %满足分配功率等式约束
          Pdzl(a)  = phen1(a,2);
          Pdc(a)   = phen1(a,1) +  phen1(a,2) - Pfd;
          
          
          %约束条件，重复约束条件
          if Pdc(a) <= A1;Pdc(a)=A1;end;
          if Pdc(a) >= B1;Pdc(a)=B1;end;
          
          if Pfl(a) <= A2;Pfl(a)=A2;end;
          if Pfl(a) >= B2;Pfl(a)=B2;end;   
          
          if Pdzl(a) <= A3;Pdzl(a)=A3;end;
          if Pdzl(a) >= B3;Pdzl(a)=B3;end; 
          
          
          
          %计算对应的目标值
          Moneys   = func_obj(Pdc(a),Pfl(a),Pdzl(a),Pfd,t,Pricefd);
          JJ(a,1)  = Moneys;
      end 
      Objvsel     = JJ ;    
      [Chrom,Objv]= reins(Chrom,Selch,1,1,Objv,Objvsel);   
      gen         = gen+1; 
      [VS,IS]     = min(Objv);
      Best(gen)   = 1/min(Objv);
      fit(gen)    = min(Objv);
      Pdc2(gen)   = Pdc(IS);
      Pfl2(gen)   = Pfl(IS);
      Pdzl2(gen)  = Pdzl(IS);
end 
     if t == 1
        figure;
        plot(Best,'b','linewidth',2);
        xlabel('迭代次数');
        ylabel('优化后的收益');
        grid on
        figure;
        plot(fit,'b','linewidth',2);
        xlabel('迭代次数');
        ylabel('适应度函数收敛情况');
        grid on
     end
     
     Timesdelay
     
    Pdc21  = [Pdc21,Pdc2(end)];
    Pfl21  = [Pfl21,Pfl2(end)];
    Pdzl21 = [Pdzl21,Pdzl2(end)];   

    if t <= Timesdelay
        Pdc22  = [Pdc22,mean(Pdc21(1:t))];
        Pfl22  = [Pfl22,mean(Pfl21(1:t))];
        Pdzl22 = [Pdzl22,mean(Pdzl21(1:t))]; 
    else
        Pdc22  = [Pdc22,mean(Pdc21(t-Timesdelay:t))];
        Pfl22  = [Pfl22,mean(Pfl21(t-Timesdelay:t))];
        Pdzl22 = [Pdzl22,mean(Pdzl21(t-Timesdelay:t))]; 
    end
end
02_043m