%风、光、负荷数据：比利时2023年2月7日
clc
clear all
close all
%风机功率波动
Wind_data=xlsread('Wind_data.xlsx');
[wind_data,ps1]=mapminmax(Wind_data',0,1);
wind_data=wind_data';
%光伏功率波动
Ph_data=xlsread('Ph_data.xlsx');
[ph_data,ps2]=mapminmax(Ph_data',0,1);
ph_data=ph_data';
%负荷波动
PQ_data=xlsread('PQ_data.xlsx');
PQ_data(97,:)=8000;
[PQ_data,ps3]=mapminmax(PQ_data',1,2);
PQ_data=PQ_data';
%PSAT文件加载
initpsat;
clpsat.readfile=0;
runpsat('IEEE33bw','data');
runpsat('pf');
%导入波动数据
for i=1:96
wind(i,1)=wind_data(i,1)*0.08;
predicted_wind(i,1)=wind_data(i,2)*0.08;
ph(i,1)=ph_data(i,1)*0.08;
predicted_ph(i,1)=ph_data(i,2)*0.08;
LD(i,:)=PQ_data(i,1)*PQ.con(1:32,4)';
predicted_LD(i,:)=PQ_data(i,2)*PQ.con(1:32,4)';
end
predicted_wind(97,1)=wind_data(1,2)*0.08;
predicted_ph(97,1)=ph_data(1,2)*0.08;
predicted_LD(97,:)=PQ_data(1,2)*PQ.con(1:32,4)';
predicted_wind(98,1)=wind_data(2,2)*0.08;
predicted_ph(98,1)=ph_data(2,2)*0.08;
predicted_LD(98,:)=PQ_data(2,2)*PQ.con(1:32,4)';
%日内无控制潮流计算
for i=1:96
PQgen.store(1:3,4)=wind(i,1);
PQgen.store(2,4)=1.5.*wind(i,1);
PQgen.store(4:5,4)=ph(i,1);
PQgen.store(5,4)=1.5.*ph(i,1);
PQ.store(:,4)=LD(i,:)';
runpsat('pf');
voltages(:,i) = DAE.y(1+Bus.n:2*Bus.n);
Uwind(i,1:3)=voltages([7;16;21],i);
Uph(i,1:2)=voltages([24;30],i);
Us(i,:)=voltages([1:6,8:15,17:20,22:23,25:29,31:33],i);
end
%绘图
plt(wind,ph,LD,Us,Uph,Uwind)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%控制%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=waitbar(1/96,'求解进度');
for i=1:96
%电网上传断面数据
PQgen.store(1:3,4)=wind(i,1);
PQgen.store(4:5,4)=ph(i,1);
PQ.store(:,4)=LD(i,:)';
runpsat('pf');
sw=SW.con;pqgen=PQgen.con;pq_store=PQ.con;
sw_store(:,:,i)=sw;pqgen_store(:,:,i)=pqgen;
voltages = DAE.y(1+Bus.n:2*Bus.n);
voltages_uncontrol(:,i)=voltages;
%扰动法灵敏度计算
[sensitivity_H_U,sensitivity_H_Ploss,...
sensitivity_Wind_U_P,sensitivity_Wind_Ploss_P,sensitivity_Wind_U_Q,sensitivity_Wind_Ploss_Q,...
sensitivity_Ph_U_P,sensitivity_Ph_Ploss_P,sensitivity_Ph_U_Q,sensitivity_Ph_Ploss_Q]=sensitivity(i,LD,pqgen,sw);
%%%%%%%%%%%%%%%%%%%%%%控制%%%%%%%%%%%%
vref=1.0;
x=sdpvar(5,1);
%单时间断面控制
v=voltages(2:33,1)+sensitivity_Wind_U_Q(:,2:33)'*x(1:3,1)+sensitivity_Ph_U_Q(:,2:33)'*x(4:5,1);
%目标函数
c1=100;c2=10;
J=c1.*sum((v-vref).^2)+c2.*sum(x);
F=[
(pqgen(:,5)+x)<=(0.15.^2-pqgen(:,4).^2).^0.5;
0.95<=v&v<=1.05;
];
options = sdpsettings('solver','cplex');
sol = optimize(F,J,options);
x_ANS=double(x);
if sol.problem==1
disp('求解出错')
break
end
PQgen.store(:,5)=pqgen(:,5)+x_ANS;
runpsat('pf');
voltages_control(:,i) = DAE.y(1+Bus.n:2*Bus.n);
str=['求解进度...',num2str(i/96*96),'/96'];
waitbar(i/96,h,str)
end
delete(h);
figure
plot(voltages_uncontrol');
hold on
plot([0,96],[1.05,1.05],'r')
hold on
plot([0,96],[0.95,0.95],'r')
legend('未控制前');
figure
plot(voltages_control');
hold on
plot([0,96],[1.05,1.05],'r')
hold on
plot([0,96],[0.95,0.95],'r')
legend('控制后');在此基础上使负荷固定比例持续增长，观察电压网损变化，避免维度错误
EstimatedSensitivity：快速灵敏度估算在FTM模式下的应用
EstimatedSensitivity是一种基于IQcapture和CW输入信号的敏感度估算方法，它在FTM模式下可用，能缩短测试时间，允许在发射机不工作时进行测试，尤其适用于debug和Txdesense测试。该算法与噪声系数和最低信噪比相关，误差通常在±1dB内，前提是spur小于25dB。
扫一扫
8万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
