目录
1 主要内容
2 部分代码
3 程序结果
4 程序链接
1 主要内容
该程序完美复现《Optimal air-conditioning load control in distribution network with intermittent renewables》，中文题目（翻译）为《含可再生能源的配电网最佳空调负荷优化控制》，实现内容：为了减轻包括风力发电、环境温度变化和电力零售价格在内的随机变量的不确定性，采用滚动水平优化（RHO）策略来连续更新实时信息并进入控制窗口。此外，为了确保客户的热舒适性，引入了一种新颖的两参数热模型来更精确地计算室内温度变化，保证用户舒适度，以系统总运行成本最低为目标，采用基于MILP（混合整数线性规划）的方法来调度可中断的空调负荷，从而实现最佳调度效果。
- 该程序采用matlab+mosek(求解器)进行求解，具体mosek安装方法可参见官网，后期公众号也会更新该求解器的安装注意事项。
- 两参数模型如下：
房屋的热过程由两个部分组成。一个组成部分是房屋内部的热质量，另一个是具有显着不同的热容量的墙壁的热质量。考虑到墙壁的热容量，房屋的室内空气温度变化可能会显著不同。由于热模型的复杂性会对冷却能量的计算精度产生重大影响，因此该模型能够虽然复杂性增加，但是准确性也有效提升。
- RHO策略
RHO的主要程序如下。
1）在第一个时间步，基于日前的预测数据确定包括电价，环境温度和风能输出在内的参数。MILP模型会根据最低运营成本目标计算并生成一组参数（例如Tr和Tw）。
2）在下一个时间步，基于更新的输入参数（包括实时价格（RTP），新近预测的环境温度，更新的未来风能输出，前一步生成的数据Tr等），使用MILP模型优化成本函数，并为下一个控制窗口生成一组新的参数。
3）在每个时间步，向前移动控制窗口，重复上述过程，直到完成计划范围的最后一个时间步。
2 部分代码
%% Variables
syms Sac1; % (integer variables) operational status of AC--0=off,1=on各个房间空调开关状态，房间温度，墙体温度
syms Tr1; % temperature of the room
syms Tw1; % temperature of the wall
syms Sac2; % (integer variables) operational status of AC--0=off,1=on
syms Tr2; % temperature of the room
syms Tw2; % temperature of the wall
syms Sac3; % (integer variables) operational status of AC--0=off,1=on
syms Tr3; % temperature of the room
syms Tw3; % temperature of the wall
syms Sac4; % (integer variables) operational status of AC--0=off,1=on
syms Tr4; % temperature of the room
syms Tw4; % temperature of the wall
syms Sac5; % (integer variables) operational status of AC--0=off,1=on
syms Tr5; % temperature of the room
syms Tw5; % temperature of the wall
syms P_bat; % power output from battery to consumers (kW)消费者从电池使用的电量
syms P_buy; % power bought from grid (kW)从电网购买的电量
syms P_sell; % power from battery to grid (kW)电池输入电网的电量
syms SoC; % battery state of charge电池电量，百分比
numofvars = 19; % the number of variables which have been listed above
%% Import ambient temperature and electrcity prices data from spreadsheets
Tamb = xlsread('Main.xlsx','Data','P2:P97')'; % importing data for expected consumption during the day (pu)预期环境温度
c_buy = xlsread('Main.xlsx','Data','D2:D97')'; % importing data for the cost of purchasing power from the grid (pu)购入电价
PCritical = xlsread('Main.xlsx','Data','Q2:Q49'); % importing data for the critical loads基础负荷，不可控部分
PCritical=0.2*Process_Load(PCritical);
WindGeneration=xlsread('Main.xlsx','Data','U1:U96')';% importing data for wind generation风力发电预测数据
WindGeneration=0.1*WindGeneration;
c_sell = xlsread('Main.xlsx','Data','E2:E97')'; % importing data for the price of selling power to the grid (kw)卖出电价
%% Data
t_step = 15; % timestep in minutes每次滑动距离（分钟）
t_total = 24*60; % total time window in minutes总时长
N = t_total / t_step; % total number of time steps总滑动次数
E_bat = 60; % battery capacity 1000 kWH电池容量
P_max_discharge = 60; % battery maximum discharge output kw 电池最大输出功率
P_max_bound = 10000;
density_of_wall=2400; % Density of wall is 2000kg/m3墙体密度
density_of_air=1.1839; % Density of air at 25℃空气密度
GroupOneNumber=10; % House number of group one is 100;每组房子的房屋数量
GroupTwoNumber=12; % House number of group two is 100;
GroupThreeNumber=14; % House number of group three is 100;
GroupFourNumber=18; % House number of group four is 100;
GroupFiveNumber=15; % House number of group five is 100;
%House 1
Tr_init_1 = 25; % initial room temperature is 25℃房间初始温度
Tamb_init_1 = 25; % last time step ambient temperature before optimization is 25℃
Tw_init_1=25; % initial wall temperature is 25℃初始墙体温度
Sac_init_1 =0; % initial ac operational status is off初始空调状态为关闭
cpa_1 = 1005.4; % cp = cp of air (273 K) = 1005.4 J/kg-K空气比热容
cpw_1=840; % cpw = cp of wall (273k)=840J/kg-k墙体比热容
Pac_1 =3; % power of air conditioner(kW)空调的功率
COP_1=2.8; % coefficient of performance空调性能系数
Qac_1=COP_1*Pac_1*1000*0.25*900; % energy of air conditioner at 15minutes(J)15分钟内空调消耗的电量
NumberOfHouse1=200; % Number of House 1 is 200
%House 2
Tr_init_2 = 25; % initial room temperature is 25℃
Tamb_init_2 = 25; % last time step ambient temperature before optimization is 25℃
Tw_init_2=25; % initial wall temperature is 25℃
Sac_init_2 =0; % initial ac operational status is off
cpa_2 = 1005.4; % cp = cp of air (273 K) = 1005.4 J/kg-K
cpw_2=840; % cpw = cp of wall (273k)=840J/kg-k
Pac_2 = 4; % power of air conditioner(kW)
COP_2= 2.8; % coefficient of performance
Qac_2=COP_2*Pac_2*1000*0.25*900; % energy of air conditioner at 15minutes(J)
%House 3
Tr_init_3 = 25; % initial room temperature is 25℃
Tamb_init_3 = 25; % last time step ambient temperature before optimization is 25℃
Tw_init_3=25; % initial wall temperature is 25℃
Sac_init_3 =0; % initial ac operational status is off
cpa_3 = 1005.4; % cp = cp of air (273 K) = 1005.4 J/kg-K
cpw_3=840; % cpw = cp of wall (273k)=840J/kg-k
Pac_3 =4; % power of air conditioner(kW)
COP_3=2.8; % coefficient of performance
Qac_3=COP_3*Pac_3*1000*0.25*900; % energy of air conditioner at 15minutes(J)
%House 4
Tr_init_4 = 25; % initial room temperature is 25℃
Tamb_init_4 = 25; % last time step ambient temperature before optimization is 25℃
Tw_init_4=25; % initial wall temperature is 25℃
Sac_init_4 =0; % initial ac operational status is off
cpa_4 = 1005.4; % cp = cp of air (273 K) = 1005.4 J/kg-K
cpw_4=840; % cpw = cp of wall (273k)=840J/kg-k
Pac_4 = 2.8; % power of air conditioner(kW)
COP_4=2.8; % coefficient of performance
Qac_4=COP_4*Pac_4*1000*0.25*900; % energy of air conditioner at 15minutes(J)
%House 5
Tr_init_5 = 25; % initial room temperature is 25℃
Tamb_init_5 = 25; % last time step ambient temperature before optimization is 25℃
Tw_init_5=25; % initial wall temperature is 25℃
Sac_init_5 =0; % initial ac operational status is off
cpa_5 = 1005.4; % cp = cp of air (273 K) = 1005.4 J/kg-K
cpw_5=840; % cpw = cp of wall (273k)=840J/kg-k
Pac_5 =3.5; % power of air conditioner(kW)
COP_5=2.8; % coefficient of performance
Qac_5=COP_5*Pac_5*1000*0.25*900; % energy of air conditioner at 15minutes(J)
%% Define the house geometry各房间的几何参数
% converst radians to degrees
r2d = 180/pi;
% Roof pitch = 40 deg屋顶坡度
pitRoof = 40/r2d;
%--House one--%
% House length = 30 m长度
lenHouse1 = 20;
% House width = 10 m宽度
widHouse1 = 12;
% House height = 4 m高度
LWall1 = 0.24;
% Width of the wall =0.24m墙体厚度
htHouse1 = 3.2;
% Number of windows = 6窗口数量
numWindows1 = 6;
% Height of windows = 1 m窗口高度
htWindows1 = 1;
% Width of windows = 1 m窗口宽度
widWindows1 = 1;
windowArea1 = numWindows1*htWindows1*widWindows1;%窗户总面积
wallArea1 = 2*lenHouse1*htHouse1 + 2*widHouse1*htHouse1 + ...
2*(1/cos(pitRoof/2))*widHouse1*lenHouse1 + ...
tan(pitRoof)*widHouse1 - windowArea1;%墙体总面积
volume_of_wall1=lenHouse1*htHouse1*LWall1*2+(widHouse1-LWall1*2)*htHouse1*LWall1*2-windowArea1*LWall1;%墙体体积
volume_of_air1=lenHouse1*htHouse1*widHouse1+0.5*widHouse1*tan(40)*0.5*widHouse1;%空气体积
Mw1=volume_of_wall1*density_of_wall;%墙体总质量
Mair1 = density_of_air*volume_of_air1; % mass of air within the household(kg)空气总质量
% Define the type of insulation used定义绝热材料
% Glass wool in the walls, 0.2 m thick墙面玻璃棉，0.2米厚度
% k is in units of J/sec/m/C - convert to J/hr/m/C multiplying by 3600
kWall1 = 0.038*900; % hour is the time unit
RWall1 = LWall1/(kWall1*wallArea1);
% Glass windows, 0.01 m thick玻璃窗，0.01米厚度
kWindow1 = 0.78*900; % hour is the time unit
LWindow1 = .01;
RWindow1 = LWindow1/(kWindow1*windowArea1);
% Determine the equivalent thermal resistance for the whole building
%确定整个建筑的热阻抗
Rw1 = RWall1*RWindow1/(RWall1 + RWindow1);
% Convection Heat transfer coefficient between air and wall,hwr=15w/(m2k);
%空气与墙壁对流换热系数hwr
hwr1=15;
% convection resistance is Rwr=1/(hA)对流阻抗Rwr
Rwr1=1/(hwr1*wallArea1);
Rwa1=Rwr1;
%Calculate the constants to simplified the计算常数以简化计算
3 程序结果
4 程序链接
点击直达！
