matlab计算最优潮流

%第一步先由内点法求得局部最优解
%%%对比用
%%内点法加SDP的全局最优性条件考虑了提前判断起作用约束的条件
clc;clear all;
mpc=loadcase（‘case9_c1‘）; %可加载不同系统
define_constants; %定义常量
baseMVA=mpc.baseMVA; %系统基值
n=size（mpc.bus1）;  %节点个数
m=size（mpc.gen1）;  %发电机个数
branchnum=size（mpc.branch1）;   %支路数

% format long
f=0; %内点法求得的发电成本
Pg=zeros（m1）; %内点法求得的机组出力
bus_slack=find（mpc.bus（:BUS_TYPE）==3）; %平衡节点
v=1; %内点法求得的平衡节点电压

    result=runopf（mpc）; %求解OPF
    
%     [baseMVAbusgengencostbranchfsuccesset]=runopf（casename）;
    
    Pg（:1）=result.gen（:PG） %机组出力
    v=result.bus（bus_slackVM）; %平衡节点电压
    f=f+result.f;
    orig_obj=0;
for k=1:m
    orig_obj=orig_obj+mpc.gencost（k5）*Pg（k）^2+mpc.gencost（k6）*Pg（k）+mpc.gencost（k7）;
end

f   %目标函数值
v   %平衡节点电压幅值
bus_ma=result.bus（:VM）;
bus_an=（result.bus（:VA）./180）*pi;

X=zeros（2*n1）;
V_V=zeros（n1）;%v-v是表示电压实部虚部写在一起的值
for k=1:n
    X（k1）=bus_ma（k）*cos（bus_an（k））;
    X（n+k1）=bus_ma（k）*sin（bus_an（k））;
    V_V（k）=（X（k1）+X（n+k1）*i）*10;
end
V_V
X=X.*10;   %注意在方程中所使用的值是标幺值还是绝对值
W=zeros（2*n2*n）;
W=X*（X.‘）;

%% 以下为验证最优解是否为全局最优解
define_constants; %定义常量

%% 生成节点导纳矩阵
% Y
Y=zeros（nn）;
for k=1:branchnum
    Y（mpc.branch（k1）mpc.branch（k2））=-1/（mpc.branch（kBR_R）+mpc.branch（kBR_X）*1i）;
    Y（mpc.branch（k2）mpc.branch（k1））=-1/（mpc.branch（kBR_R）+mpc.branch（kBR_X）*1i）;
    C（mpc.branch（k1）mpc.branch（k2））=mpc.branch（kBR_B）*1i/2; %charging capacity
    C（mpc.branch（k2）mpc.branch（k1））=mpc.branch（kBR_B）*1i/2;
end
for k=1:n
    Y（kk）=-sum（Y（k:））+sum（C（k:））;
end
for k=1:branchnum
    if mpc.branch（kTAP）~=0%turns ratio
        x2=mpc.branch（kF_BUS）;
        x1=mpc.branch（kT_BUS）;
        a_aux=mpc.branch（kTAP）;
        Y（x1x2）=Y（x1x2）+1/（mpc.branch（k4）*1i）-1/（mpc.branch（k4）*1i）/a_aux;
        Y（x2x2）=Y（x2x2）-1/（mpc.branch（k4）*1i）+1/（mpc.branch（k4）*1i）/a_aux^2;
        Y（x2x1）=Y（x2x1）+1/（mpc.branch（k4）*1i）-1/（mpc.branch（k4）*1i）/a_aux;
    end
end

I=Y*V_V;%每个节点的注入电流
Power=V_V.*transpose（I‘）;%每个节点的注入功率，电流需要共轭转置
Power=[Power Power];
Power（:1）=real（Power（:1））+mpc.bus（:PD）;
Power（:2）=imag（Power（:2））+mpc.bus（:QD）;
% [[1:n]‘ Power] %各节点的发电机出力，无发电机的节点值为0
Power





V_V;
V_rad=zeros（2*n1）;
for k=1:n
    V_rad（k1）=abs（V_V（k））;
    V_rad（n+k1）=angle（V_V（k））;
end
Power;
X_rad=zeros（2*m+2*n1）;
for k=1:m
    X_rad（k1）=real（Power（mpc.gen（k1）1））;
    X_rad（m+k1）=（Power（mpc.gen（k1）2））;
end
for k=1:n
    X_rad（2*m+k1）=abs（V_V（k））;
    X_rad（2*m+n+k1）=angle（V_V（k））;
end
% 下面为变量初始化
s_0=zeros（2*m+n+branchnum1）;
z_0=zeros（2*m+n+branchnum1）;
lambuda_0=zeros（2*n1）;
pai_0=zeros（2*m+n+branchnum1）;
v_0=zeros（2*m+n+branchnum1）;

gama_s=zeros（2*m+n+branchnum1）;
gama_z=zeros（2*m+n+branchnum1）;
gama_pai=zeros（2*m+n+branchnum1）;
gama_v=zeros（2*m+n+branchnum1）;
gama_x=zeros（2*m+2*n1）;
gama_la=zeros（2*n1）;

gama=0;
mune0=0.1;
for k=1:m
    s_0（k）=min（max（gama*（mpc.gen（k9）-mpc.gen（k10））X_rad（k）-mpc.gen（k10））