【优化电价】基于matlab内点法求解实时电价最优问题【含Matlab源码 1161期】
        【摘要】 
                    
                        
                    
                    一、内点法简介 
 法。而且无论是面对LP还是QP,内点法都显示出了相当的极好的性能,例如多项式的算法复杂度。 
二、部分源代码 
clear;
clc;
errArr=[];
%%
%初始化!!!
in...
    
    
    
    一、内点法简介

 法。而且无论是面对LP还是QP,内点法都显示出了相当的极好的性能,例如多项式的算法复杂度。
二、部分源代码
clear;
clc;
errArr=[];
%%
%初始化!!!
initial;
% Start clock
t1 = clock;
%%
ROU=sl'*MU_MIN+su'*MU_MAX;
MUt=SIGMA*ROU/(2*length(sl));%初始对偶因子与惩罚因子计算%
ik=0;%计迭代次数!!!
%迭代循环过程!!
while(abs(ROU)>=err)
    %%
    %Calcute h,g matrix
    ROU=sl'*MU_MIN+su'*MU_MAX;
        errArr=[errArr;ROU;];
    SIGMA=0;
    MU=SIGMA*ROU/(2*length(sl));         %中心参数置零%
    for i=1:30
        temp=0;
        for j=1:30
            temp=temp-V(j)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));
        end
        if (i>6) 
            tPg=0;
        else
            tPg=Pg(i);
        end
        h(i)=tPg-Pd(i)+V(i)*temp;
    end
    for i=1:30
        temp=0;
        for j=1:30
            temp=temp-V(j)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));
        end
        if (i>6) 
            tQg=0;
        else
            tQg=Qg(i);
        end
        h(i+30)=tQg-Qd(i)+V(i)*temp;
    end        % Cal  h END
    
    for i=1:6
        g(i)=Pg(i);
        g(i+6)=Qg(i);
    end
    for i=1:30
        g(i+12)=V(i);
    end        % Cal  g  END
    %Calcute h,g matrix END
    %%
    %Calculate Jacobian&Hessian matix
    %First Step: Jf,Hf
    for i=1:6
        Jf(i)=2*gencost(i,5)*Pg(i)+gencost(i,6);
        Hf(i,i)=2*gencost(i,6);
    end
    %Second Step: Jh, h为等式约束
    for i=1:6         %前6行对Pg求导,由此已求出
        Jh(i,i)=1; 
    end
    for i=7:12        %7-12行对Qg求导,由此已求出
        Jh(i,i+24)=1; 
    end
    for i=1:30       %形成13-42行的1-60列
        for j=1:30
            tempVp=0;
            tempVq=0;
            if (j==i)
                for k=1:30
                    tempVp=tempVp-V(k)*aY(j,k)*cos(Vth(j)-Vth(k)-Yth(j,k));
                    tempVq=tempVq-V(k)*aY(j,k)*sin(Vth(j)-Vth(k)-Yth(j,k));
                end
                Jh(12+j,i)=tempVp-aY(j,j)*V(j)*cos(Yth(j,j));
                Jh(12+j,30+i)=tempVq+aY(j,j)*V(j)*sin(Yth(j,j));
            else
                Jh(12+j,i)=-aY(i,j)*V(i)*cos(Vth(i)-Vth(j)-Yth(i,j));
                Jh(12+j,30+i)=-aY(i,j)*V(i)*sin(Vth(i)-Vth(j)-Yth(i,j));
            end
        end
    end
    for i=1:30       %形成43-72行的1-60列
        for j=1:30
            tempVp=0;
            tempVq=0;
            if (j==i)
                for k=1:30
                    tempVp=tempVp+aY(j,k)*V(k)*sin(Vth(j)-Vth(k)-Yth(j,k));
                    tempVq=tempVq-aY(j,k)*V(k)*cos(Vth(j)-Vth(k)-Yth(j,k));
                end
                tempVp=tempVp-V(j)*aY(j,j)*sin(-Yth(j,j));
                tempVq=tempVq+V(j)*aY(j,j)*cos(-Yth(j,j));
                Jh(42+j,i)=V(i)*tempVp;
                Jh(42+j,30+i)=V(i)*tempVq;
            else
                Jh(42+j,i)=-aY(i,j)*V(i)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j));
                Jh(42+j,30+i)=aY(i,j)*V(i)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j));
            end
        end
    end
    %Third Step: Hh
    %有功部分
    for i=1:30
        for j=1:30
            for k=j:30
                if (j==k&&i~=j) 
                    Hh(j+12,k+12,i)=0;          %VV
                    Hh(j+42,k+42,i)=V(i)*aY(i,j)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j)); %%thth
                elseif (j==k&&i==j)
                    Hh(j+12,k+12,i)=-2*aY(j,j)*cos(Yth(i,i));  %VV
                    temp=0;                   %thth
                    for l=1:30
                        temp=temp+aY(j,l)*V(l)*cos(Vth(j)-Vth(l)-Yth(j,l));
                    end
                    temp=temp-aY(i,i)*V(i)*cos(-Yth(i,i));
                    Hh(j+42,k+42,i)=V(i)*temp;
                elseif (k==i)
                    Hh(j+12,k+12,i)=-aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));    %VV
                    Hh(k+12,j+12,i)=Hh(j+12,k+12,i);
                    Hh(j+42,k+42,i)=-V(i)*aY(i,j)*V(j)*cos(Vth(i)-Vth(j)-Yth(i,j));  %thth
                    Hh(k+42,j+42,i)=Hh(j+42,k+42,i);
                elseif (j==i)
                    Hh(j+12,k+12,i)=-aY(i,k)*cos(Vth(i)-Vth(k)-Yth(i,k));    %VV
                    Hh(k+12,j+12,i)=Hh(j+12,k+12,i);
                    Hh(j+42,k+42,i)=-V(i)*aY(i,k)*V(k)*cos(Vth(i)-Vth(k)-Yth(i,k));  %thth
                    Hh(k+42,j+42,i)=Hh(j+42,k+42,i);
                end
            end
        end
    end                 %至此已形成(13-42,13-42)和(42-72,43-72)
    for i=1:30
        for j=1:30
            for k=1:30
                if (j==k&&i~=j) 
                    Hh(j+42,k+12,i)=-V(i)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));      %thV
                elseif (j==k&&i==j)
                    temp=0;                   %thV
                    for l=1:30
                        temp=temp+aY(j,l)*V(l)*sin(Vth(j)-Vth(l)-Yth(j,l));    
                    end
                    Hh(j+42,k+12,i)=temp-V(i)*aY(i,i)*sin(-Yth(i,i));
                elseif (j==i)
                    Hh(j+42,k+12,i)=V(i)*aY(i,k)*sin(Vth(i)-Vth(k)-Yth(i,k));   %thV   
                elseif (k==i)
                    Hh(j+42,k+12,i)=-V(j)*aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));  %thV
                end
            end
        end
        Hh(13:42,43:72,i)=Hh(43:72,13:42,i)';
    end                           %至此已形成(42-72,13-42)和(13-42,43-72)
    %无功部分
    for i=1:30
        for j=1:30
            for k=j:30
                if (j==k&&i~=j) 
                    Hh(j+12,k+12,i+30)=0;          %VV
                    Hh(j+42,k+42,i+30)=V(i)*aY(i,j)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j)); %%thth
                elseif (j==k&&i==j)
                    Hh(j+12,k+12,i+30)=2*aY(j,j)*sin(Yth(i,i));  %VV
                    temp=0;                   %thth
                    for l=1:30
                        temp=temp+aY(j,l)*V(l)*sin(Vth(j)-Vth(l)-Yth(j,l));
                    end
                    temp=temp-aY(i,i)*V(i)*sin(-Yth(i,i));
                    Hh(j+42,k+42,i+30)=V(i)*temp;
                elseif (k==i)
                    Hh(j+12,k+12,i+30)=-aY(i,j)*sin(Vth(i)-Vth(j)-Yth(i,j));    %VV
                    Hh(k+12,j+12,i+30)=Hh(j+12,k+12,i+30);
                    Hh(j+42,k+42,i)=-V(i)*aY(i,j)*V(j)*sin(Vth(i)-Vth(j)-Yth(i,j));  %thth
                    Hh(k+42,j+42,i+30)=Hh(j+42,k+42,i+30);
                elseif (j==i)
                    Hh(j+12,k+12,i+30)=-aY(i,k)*sin(Vth(i)-Vth(k)-Yth(i,k));    %VV
                    Hh(k+12,j+12,i+30)=Hh(j+12,k+12,i+30);
                    Hh(j+42,k+42,i+30)=-V(i)*aY(i,k)*V(k)*sin(Vth(i)-Vth(k)-Yth(i,k));  %thth
                    Hh(k+42,j+42,i+30)=Hh(j+42,k+42,i+30);
                end
            end
        end
    end                 %至此已形成(13-42,13-42)和(42-72,43-72)
    for i=1:30
        for j=1:30
            for k=1:30
                if (j==k&&i~=j) 
                    Hh(j+42,k+12,i+30)=V(i)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));      %thV
                elseif (j==k&&i==j)
                    temp=0;                   %thV
                    for l=1:30
                        temp=temp-aY(j,l)*V(l)*cos(Vth(j)-Vth(l)-Yth(j,l));    
                    end
                    Hh(j+42,k+12,i+30)=temp+V(i)*aY(i,i)*cos(-Yth(i,i));
                elseif (j==i)
                    Hh(j+42,k+12,i+30)=-V(i)*aY(i,k)*cos(Vth(i)-Vth(k)-Yth(i,k));   %thV   
                elseif (k==i)
                    Hh(j+42,k+12,i+30)=V(j)*aY(i,j)*cos(Vth(i)-Vth(j)-Yth(i,j));  %thV
                end
            end
        end
        Hh(13:42,43:72,i+30)=Hh(43:72,13:42,i+30)';
    end                           %至此已形成(42-72,13-42)和(13-42,43-72)
    %Hh形成完毕
    %Fourth Step: Jg, Hg
    Jg=eye(42,42);
    Jg=[Jg;zeros(30,42)];
    Hg=zeros(72);
    %Calculation Jacobian&Hessian matrix END
    %%
    %Calculate Newton Iteration 误差迭代量
    %Cal LX0-------------------------1
    LX0=Jf-Jh*Lam+Jg*(-MU_MIN+MU_MAX);
    %Cal LLam-------------------------2
    LLam0=h;
    pferr=max(LLam0);
    %Cal LMU_MIN-------------------------3
    LMU_MIN0=g-sl-gmin;
    %Cal LMU_MAX-------------------------4
  
 - 1
 - 2
 - 3
 - 4
 - 5
 - 6
 - 7
 - 8
 - 9
 - 10
 - 11
 - 12
 - 13
 - 14
 - 15
 - 16
 - 17
 - 18
 - 19
 - 20
 - 21
 - 22
 - 23
 - 24
 - 25
 - 26
 - 27
 - 28
 - 29
 - 30
 - 31
 - 32
 - 33
 - 34
 - 35
 - 36
 - 37
 - 38
 - 39
 - 40
 - 41
 - 42
 - 43
 - 44
 - 45
 - 46
 - 47
 - 48
 - 49
 - 50
 - 51
 - 52
 - 53
 - 54
 - 55
 - 56
 - 57
 - 58
 - 59
 - 60
 - 61
 - 62
 - 63
 - 64
 - 65
 - 66
 - 67
 - 68
 - 69
 - 70
 - 71
 - 72
 - 73
 - 74
 - 75
 - 76
 - 77
 - 78
 - 79
 - 80
 - 81
 - 82
 - 83
 - 84
 - 85
 - 86
 - 87
 - 88
 - 89
 - 90
 - 91
 - 92
 - 93
 - 94
 - 95
 - 96
 - 97
 - 98
 - 99
 - 100
 - 101
 - 102
 - 103
 - 104
 - 105
 - 106
 - 107
 - 108
 - 109
 - 110
 - 111
 - 112
 - 113
 - 114
 - 115
 - 116
 - 117
 - 118
 - 119
 - 120
 - 121
 - 122
 - 123
 - 124
 - 125
 - 126
 - 127
 - 128
 - 129
 - 130
 - 131
 - 132
 - 133
 - 134
 - 135
 - 136
 - 137
 - 138
 - 139
 - 140
 - 141
 - 142
 - 143
 - 144
 - 145
 - 146
 - 147
 - 148
 - 149
 - 150
 - 151
 - 152
 - 153
 - 154
 - 155
 - 156
 - 157
 - 158
 - 159
 - 160
 - 161
 - 162
 - 163
 - 164
 - 165
 - 166
 - 167
 - 168
 - 169
 - 170
 - 171
 - 172
 - 173
 - 174
 - 175
 - 176
 - 177
 - 178
 - 179
 - 180
 - 181
 - 182
 - 183
 - 184
 - 185
 - 186
 - 187
 - 188
 - 189
 - 190
 - 191
 - 192
 - 193
 - 194
 - 195
 - 196
 - 197
 - 198
 - 199
 - 200
 - 201
 - 202
 - 203
 - 204
 - 205
 - 206
 - 207
 - 208
 - 209
 - 210
 - 211
 - 212
 - 213
 - 214
 - 215
 - 216
 - 217
 - 218
 
三、运行结果

 
 
 
四、matlab版本及参考文献
1 matlab版本
 2014a
2 参考文献
 [1] 包子阳,余继周,杨杉.智能优化算法及其MATLAB实例(第2版)[M].电子工业出版社,2016.
 [2]张岩,吴水根.MATLAB优化算法源代码[M].清华大学出版社,2017.
文章来源: qq912100926.blog.csdn.net,作者:海神之光,版权归原作者所有,如需转载,请联系作者。
原文链接:qq912100926.blog.csdn.net/article/details/119220848
        【版权声明】本文为华为云社区用户转载文章,如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱:
            cloudbbs@huaweicloud.com
        
        
        
        
        
        
        - 点赞
 - 收藏
 - 关注作者
 
            
           
评论(0)