%第一步先由內點法求得局部最優解
%%%對比用
%%內點法加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））