%第一步先由內(nèi)點(diǎn)法求得局部最優(yōu)解
%%%對(duì)比用
%%內(nèi)點(diǎn)法加SDP的全局最優(yōu)性條件考慮了提前判斷起作用約束的條件
clc;clear?all;
mpc=loadcase（‘case9_c1‘）;?%可加載不同系統(tǒng)
define_constants;?%定義常量
baseMVA=mpc.baseMVA;?%系統(tǒng)基值
n=size（mpc.bus1）;??%節(jié)點(diǎn)個(gè)數(shù)
m=size（mpc.gen1）;??%發(fā)電機(jī)個(gè)數(shù)
branchnum=size（mpc.branch1）;???%支路數(shù)

%?format?long
f=0;?%內(nèi)點(diǎn)法求得的發(fā)電成本
Pg=zeros（m1）;?%內(nèi)點(diǎn)法求得的機(jī)組出力
bus_slack=find（mpc.bus（:BUS_TYPE）==3）;?%平衡節(jié)點(diǎn)
v=1;?%內(nèi)點(diǎn)法求得的平衡節(jié)點(diǎn)電壓

????result=runopf（mpc）;?%求解OPF
????
%?????[baseMVAbusgengencostbranchfsuccesset]=runopf（casename）;
????
????Pg（:1）=result.gen（:PG）?%機(jī)組出力
????v=result.bus（bus_slackVM）;?%平衡節(jié)點(diǎn)電壓
????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???%目標(biāo)函數(shù)值
v???%平衡節(jié)點(diǎn)電壓幅值
bus_ma=result.bus（:VM）;
bus_an=（result.bus（:VA）./180）*pi;

X=zeros（2*n1）;
V_V=zeros（n1）;%v-v是表示電壓實(shí)部虛部寫(xiě)在一起的值
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;???%注意在方程中所使用的值是標(biāo)幺值還是絕對(duì)值
W=zeros（2*n2*n）;
W=X*（X.‘）;

%%?以下為驗(yàn)證最優(yōu)解是否為全局最優(yōu)解
define_constants;?%定義常量

%%?生成節(jié)點(diǎn)導(dǎo)納矩陣
%?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;%每個(gè)節(jié)點(diǎn)的注入電流
Power=V_V.*transpose（I‘）;%每個(gè)節(jié)點(diǎn)的注入功率，電流需要共軛轉(zhuǎn)置
Power=[Power?Power];
Power（:1）=real（Power（:1））+mpc.bus（:PD）;
Power（:2）=imag（Power（:2））+mpc.bus（:QD）;
%?[[1:n]‘?Power]?%各節(jié)點(diǎn)的發(fā)電機(jī)出力，無(wú)發(fā)電機(jī)的節(jié)點(diǎn)值為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））