基于PSO粒子群算法的MPPT最大功率跟踪Simulink仿真,PSO采用S函数实现

1.算法描述 MPPT控制器的全称是“最大功率点跟踪”（Maximum Power Point Tracking）太阳能控制器，是传统太阳能充放电控制器的升级换代产品。MPPT控制器能够实时侦测太阳能板的发电电压，并追踪最高电压电流值(VI)，使系统以最大功率输出对蓄电池充电。应用于太阳能光伏系统中，协调太阳能电池板、蓄电池、负载的工作，是光伏系统的大脑。

   最大功率点跟踪系统是一种通过调节电气模块的工作状态，使光伏板能够输出更多电能的电气系统能够将太阳能电池板发出的直流电有效地贮存在蓄电池中，可有效地解决常规电网不能覆盖的偏远地区及旅游地区的生活和工业用电，不产生环境污染。

    目前，光伏阵列的最大功率点跟踪(MPPT)技术，国内外已有了一定的研究，发展出各种控制方法常，常用的有一下几种：恒电压跟踪法(ConstantVoltageTracking简称CVT)、干扰观察法(PerturbationAndObservationmethod简称P&O)、增量电导法(IncrementalConductancemethod简称INC)、基于梯度变步长的电导增量法等等。(这些算法只能用在无遮挡的条件下)

1）单峰值功率输出的MPPT的算法

   目前，在无遮挡条件下，光伏阵列的最大功率点跟踪(MPPT)的控制方法常用的有以下几种：

l恒电压跟踪法(ConstantVoltageTracking简称CVT) l干扰观察法(PerturbationAndObservationmethod简称P&O) l增量电导法(IncrementalConductancemethod简称INC) l基于梯度变步长的电导增量法，等等。

2）多峰值功率输出MPPT算法

   普通的最大功率跟踪算法，如扰动观测发和电导增量法在一片云彩的遮挡下就有可能失效，不能实现真正意义的最大功率跟踪。目前，国际上也有人提出了多峰值的MPPT算法，主要包含如下三种：

结合常规算法的复合MPPT算法 Fibonacci法 短路电流脉冲法

PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己的速度和位置。

第①部分称为【记忆项】，表示上次速度大小和方向的影响；

第②部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；

第③部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。

以上面两个公式为基础，再来看一个公式：

2.仿真效果预览 matlab2022a仿真如下：

3.MATLAB核心程序

.............................
 
 
persistent Pbest;
%persistent Pbestval;
 
persistent best_index;
 
persistent c1;
persistent c2;
persistent r1;
persistent r2;
 
switch flag,
    case 0,
        sizes = simsizes;
        sizes.NumContStates = 0;
        sizes.NumDiscStates = 0;
        sizes.NumOutputs     = 1;
        sizes.NumInputs      = 2;
        sizes.DirFeedthrough = 0;
        sizes.NumSampleTimes = 1;
        sys = simsizes(sizes);
        x0=[];
        str=[];
        ts=[0.004 ,0.001];
 
        % initialize the static variables
        first = 1;
        stop = 0;
        i = 0;
        mg = 0;
        count = 3;
        Nc_count = 0;
        Nre_count = 0;
        D_out = zeros(1, NP);                       %初始化D_out、U、fitval_current、fitval_new为全0
        D_out_current = zeros(1,NP);
        best_index = 1;
        %U = zeros(1, NP);
        %V = zeros(1, NP);
        fitval_current = zeros(1, NP);
        fitval_new = zeros(1, NP);
        
        fit_order = zeros(1, NP);
        
 
        Pbest = zeros(1, NP);
        Gbest = 0;
        %Gmaxval = 0;
        %X(:) = unifrnd(XL, XU, 1,NP)
        D_out_current(:) = linspace(XL+0.005,XU-0.005,NP);                 %在[-1,1]间均匀取值
        D_out(:) = D_out_current(:);
        MoveStep = unifrnd(-MaxStep, MaxStep, 1,NP);   %初始化单个细菌的移动速度
        MoveStep_PSO = zeros(1, NP);
        stable_flag = 1;                               %初始化，电路未稳定
        stable_count = 0;
        
        c1 = 0.030;
        c2 = 0.030;
        r1 = rand();
        r2 = rand();
 
    case 3,
        %if count == 3
         %count = count + 1;
            count = 1;
        %迭代完成，输出最优值
        if stop == 1;
            [best, best_index]= max(fitval_current);
            sys = D_out(best_index);
            return;
        end          
        %判断是否迭代完成
        if mg > maxgen
            stop = 1;
        end
        
        if mg == 0  % 第一代，只进行迭代，然后计算各自功率,作为初始比较功率fitval_current   
            
            if i == 0
                i= i + 1;
                sys = D_out_current(1);  
            elseif i > 1 && i < NP
                fitval_current(i-1) = adaptfunc(D_out_current(i-1),Uin);   %记录第一代的上一个个体的功率
                Pbest(i-1) = D_out_current(i-1);
                i= i + 1;
                sys = D_out_current(i);            
            elseif i == NP
                fitval_current(i-1) = adaptfunc(D_out_current(i-1),Uin);
                Pbest(i-1) = D_out_current(i-1);
                i= i + 1;
                sys = D_out_current(NP);
            elseif i == NP + 1
                fitval_current(i-1) = adaptfunc(D_out_current(i-1),Uin);
                Pbest(i-1) = D_out_current(i-1);
                [best, best_index]= max(fitval_current);
                Gbest = D_out_current(best_index);
                 i = 1;
                 Nc_count = Nc_count + 1;
                 Nre_count = Nre_count + 1;
                 mg = mg + 1;
                 for j = 1 : NP                       %第一代种群中个体进行游动
                    MoveStep_PSO = c1*r1*(Pbest(j)-D_out_current(j)) + c2*r2*(Gbest-D_out_current(j));
                    D_out(j) = D_out_current(j)+ MoveStep(j) + MoveStep_PSO;
                    if D_out(j) <= XL
                            D_out(j) = XL+0.001;
                        elseif D_out(j) >= XU
                            D_out(j) = XU-0.001;
                    end
                 end
                 sys = D_out(1);            
            else
                i= i + 1;
                sys = D_out_current(i);            
            end
            
            return;
            
        else   %第一代之后      
            
            if i > 1 && i < NP
                fitval_new(i-1) =adaptfunc(D_out(i-1), Uin);      %计算前一次占空比下的功率
                Pbest(i-1) = D_out(i-1);
                i = i + 1;
                sys = D_out(i);
            elseif i == NP
                fitval_new(i-1) =adaptfunc(D_out(i-1), Uin);
                Pbest(i-1) = D_out(i-1);
                i = i + 1;
                sys = D_out(NP);
            elseif i == 1
                i = i + 1;            
                sys = D_out(i);       
            elseif i == NP + 1 %种群迭代一次结束，开始进行新的一代种群繁殖                            
                %趋向性操作,判定运动方向
                fitval_new(i-1) =adaptfunc(D_out(i-1), Uin);
                Pbest(i-1) = D_out(i-1);
                r1 = rand();
                r2 = rand();
                for j = 1 : NP
                if (fitval_new(j) >= fitval_current(j))
                    D_out_current(j) = D_out(j);
                    fitval_current(j) =  fitval_new(j);
                else
                    MoveStep(j) = - MoveStep(j);  %若该方向功率未改进，说明不适应生存，改变方向运动
                end       
                end
                
                c1 = c1/3;
                c2 = c2/3;
 
                 D_out(:)                              %打印最新的占空比
                [best, best_index]= max(fitval_current);  %计算最大功率点，best为最大功率，best_index为最大功率点在种群中的位置            
                mg
                D_out(best_index)                   %显示最大功率点的占空比
                %Gbest = U(best_index)               %显示最大功率点的电压
                %Gmaxval = best                      %显示最大功率点的功率
                best
                fit_order = order(fitval_current,NP)    %显示当前种群功率从大到小顺序
 
                Gbest = D_out(best_index);
                % mutation
                if Nc_count == Nc %进行满Nc次趋向性操作            
                    Nc_count = 0;
                    Nre_count = Nre_count + 1;
                    %复制操作
                    MaxStep = MaxStep/2;
                    MoveStep = unifrnd(-MaxStep,MaxStep,1,NP);                 %每Nc次繁殖生成一次新的随机步长
                    [fitval_current,D_out_current]=Reproduction(fitval_current,D_out_current,NP);      %复制
                    D_out=D_out_current;
                    if Nre_count == Nre %进行满Nre次复制操作后迁移操作
                        Nre_count = 0;
                        for j=1:NP
                        if(rand(0,1)<Ped)
                            fitval_current = 0;
                            %MoveStep(i) = unifrnd(-MaxStep,MaxStep);
                            D_out(j) = unifrnd(XL,XU);
                        end
                        end
                    end
                else
                    Nc_count = Nc_count + 1;
                    for j=1:NP
                        MoveStep_PSO = c1*r1*(Pbest(j)-D_out(j)) + c2*r2*(Gbest-D_out(j));
                        D_out(j) = D_out(j) + MoveStep(j) + MoveStep_PSO;
                        if D_out(j) <= XL
                            D_out(j) = XL+0.001;
                        elseif D_out(j) >= XU
                            D_out(j) = XU-0.001;
                        end
                    end                
                end          
            %为下一次迭代做准备
            sys = D_out(1);
            i = 1;
            mg = mg +1;
            end
        end
        return;
      
    case 2,
        Uin = u;
    case {1,4,9},
        sys = [];
end
 
A85