需求響應|動态冰蓄冷系統與需求響應政策的優化研究（Matlab代碼實作）

【現在公衆号名字改為：荔枝科研社】

📋📋📋本文目錄如下：⛳️⛳️⛳️

目錄

​​1 概述​​

​​ 2 控制算法講解​​

​​2.1 随機開環最優控制 (OLOC)​​

​​2.2 模型預測控制(MPC)​​

​​3 數學模型建立 ​​

​​3.1 實體模型​​

​​3.2 經濟模型​​

​​3.3 随機最優控制問題​​

​​4 算例及運作結果 ​​

​​5 Matlab代碼實作 ​​

​​6 寫在最後​​

1 概述

集中式空調系統作為現代城市大型建築不可缺少的部分，其耗電量越來越大，部分大中城市夏

季中央空調的耗電量已占其高峰時段用電量的 20 %以上

[1-2]

，與工業用電負荷一起造成了電力尖

峰負荷，加大了電網的峰谷負荷差，導緻白天用電 緊張、而夜晚大量電力卻無法消耗。采用冰蓄冷系統，通過在夜間低谷時段蓄冷，在白天供冷時段将冷量釋放，可以減少空調制冷機組在用電高峰期的耗能，一方面可以有效實作電力“移峰填谷”，另一方面由于存在峰谷電價差，冰蓄冷系統也節約空調整體運作費用，帶來經濟效益

。

在現階段全球氣溫升高的背景下，我覺得這個課題很值得研究。由于二氧化碳的大量排放，全球變暖。

(1)所有新型電力系統的提出思維獨到，大力發展風、光等新能源尤為重要，是以我寫了很多微電網的文章代碼。搜尋方法如下。

（2）需求響應主要分為基于價格和基于激勵的需求響應。以往電力系統排程主要靠發電側，借鑒了漂亮國，國家也開始重視需求響應。需求響應計劃鼓勵消費者在電網高峰時削減負荷。我也寫了很多關于需求響應的文章代碼。 

(3)本文的創新點：具有冰蓄冷系的建築物主要使用它來轉移冷負荷。當電價低時生産冰或冷凍水，并在電價高時儲存以給使用者提供冷能。本文用随機最優控制和模型預測控制 (MPC) 結合求解。

 2 控制算法講解

本文使用了兩種算法： 随機開環最優控制 (SOOC) 和模型預測控制 (MPC) 。

2.1 随機開環最優控制 (OLOC)

對于擾動随機的随機最優控制問題，本文采用文獻[8]中的随機最優開環控制方法(SOOC)解決。下面重點講解模型預測控制（MPC）。

2.2 模型預測控制(MPC)

2.2.1 概述

MPC 最早起源于石油化工等工業界，旨在解決經典的比例-積分-微分（PID）控制難以處理多變量限制優化控制的問題。作為一種基于模型的控制算法，其滾動優化問題中的系統模型可以是描述系統動态行為的任意形式的模型，是以從原理上該方法可處理時變或非時變、線性或非線性、有時滞或無時滞的系統限制最優控制問題。正是基于該特點，MPC 作為一類适用性非常強的控制算法在電力系統優化控制領域得到了廣泛的應用。

随着 MPC 理論日趨成熟，國内外很多企業專門研發了商業化MPC 軟體包。

模 型 預 測 控 制（model predictive control，MPC）是 近 年 來 最 引 人 關 注 的 一 類 反 饋 控 制 政策，主要由模型預測、滾動優化和回報校正 3 個環節構成。

2.2.2 作用機理

MPC 從 20 世紀七八十年代開始萌芽，經曆了模型算法控制、動态矩陣控制、廣義預測控制 3 個階段的過渡，已經發展為成熟的控制理論［21］，其原理如圖 1 所示，求解步驟為：在每一個采樣時刻，根據目前測量資訊，線上求解一個有限時域開環優化問題，并将得到的控制序列的第 1 個元素作用于被控對象。在下一個采樣時刻，重複上述過程，用新的測量值重構優化問題并重新求解［22］，可以得出如下結論。

1）與正常最優控制方法不同，MPC 的優化目标并非全局不變，而是采用滾動向前式的有限時域優化目标，以局部最優解逐級逼近全局最優解。這種有限時域優化政策可兼顧模型失配、時變等不确定性，是最優控制和不确定性情況下的折中。

2）MPC 将控制問題轉化為有限時域的滾動式開環優化問題，能考慮被控對象的未來行為特征，易于處理多耦合系統，故可結合各種優化求解算法而易于拓展。

                                                                  圖1 MPC原理

3 數學模型建立 

本文詳細對成本進行模組化。實體模型是建立在文獻[6]中提出的實體模型的基礎上，擴充到包括第二台制冰機、考慮制冰機的爬坡限制、時變性能系數以及非理想的存儲和熱交換器效率。得到的模型是一個具有線性時變動力學的完全觀測MIMO随機系統。

在後文中，變量k對集合進行了索引

={0, . . . , N− 1},N =

 / Δt為控制視界T中長度Δt (小時)的離散時間步數（本文算例中Δt=0.5）。在本文中，我們取T = 24小時，雖然稍作修改，視界可以延長到一個月或一個冷卻季.

3.1 實體模型

3.1.1 系統的狀态變量

:冰庫的充電狀态

,上界為

（本文算例中冰罐容量

為1760

）

中包含下标“th”以強調測量的是熱能，而不是電能(kWh)我們在整篇文章中都采用了這個約定。冷水機排出的熱能和所需的電能與冷水機的性能系數有關。例如，典型的制冰冷水機的性能系數為562.5-3，是以蓄熱成本為14-20$ /kWhth大緻相當于35-60 $/kWh的電力存儲成本。

:為建築冷卻不足(kWhth）即前階段所有未滿足的冷卻負荷的總和（滞後冷卻負荷）。請注意，

可能為負數，例如，如果建築物在一夜之間預冷。

3.1.2 系統的控制變量 

:配置設定給制冰的電功率。（KW）

（本文最大制冰機功率

,最大制冰機斜坡

) 

:融化冰所滿足的冷負荷。

(本文每個階段的最大融冰量

，最大融冰爬坡

)

:主冰水機組消耗的功率（KW) ，為主冷水機組（主冷卻器）消耗的功率 (kW)

(本文最大主制冰機功率為

，最大主制冰機爬坡

)

 對于每個

，i ∈{1,2,3}，定義最大工作點

和最大斜坡 

。為了強制執行斜坡限制，我們将 u(k-1) 附加到狀态，定義

。這給出了狀态限制 x(k) ∈ X =

×R4 和控制限制 u(k) ∈U(x(k))，其中:

3.1.3 系統的擾動變量 

w1(k)為建築物對冰需求( kWth )

w2 ( k )為建築物的剩餘電氣需求( kW )：照明、插頭負載等。（也可以了解為建築物的固定荷載）

我們假設擾動

是服從正态分布的。

我們還假設,對于所有 k ∈T,建築營運商在

和

的第0階段都有準确的知識。給定這些定義，狀态就按照這個下式演化：    

或者更簡潔書寫為：

其中，β∈[  0,1 ]為冰滞留率，η∈[  0,1 ]為水箱與建築物的換熱效率，

和

分别為制冷機組和主制冷機組的性能系數。時間依賴性允許冷水機組性能系數随室外溫度變化。我們假設控制器在第k階段完全了解

，盡管可以添加卡爾曼濾波器來估計噪聲或未測量狀态。 

(本文冰滞留率β取為：0.98，水箱與建築物的換熱效率）

3.2 經濟模型

3.2.1 電能成本

設 

為電價 ($/kWh)，

為淨耗電量，能源成本為：                

每天的能源價格

在前一天下午 4 點釋出，是以對于 24 小時模拟，它們是确定性已知的.

對于

，令

為受第 i 層需求費用影響的階段集，令

為第 i 層需求費用之前消耗的最大功率到階段 k，即在

∩ {0,...,k − 1} 上。變量

使用當月到目前為止的第 i 個峰值需求進行初始化，或者，如果 k = 0 是當月的開始，則使用基于曆史資料的目标需求限制。每個

服從動态:                    

a.本文取初始需求高峰(所有小時)

b.初始需求高峰(（上午 8 點至晚上 10 點）

c.初始需求高峰(（上午 8 點至晚上 6 點）

3.2.2 價格型需求響應成本

一段時間内客戶用電量波動率與價格波動率的比值稱為電價伸縮系數:

如果實施峰平谷分時電價政策，由于電力使用者對電價的反應，不同時期的用電量将發生轉移。在本文中，本文采用需求側響應電價伸縮系數矩陣描述電力需求的變化。

需求側響應電價伸縮系數矩陣由電能變化與價格變化之比的彈性系數組成。使用者用電量不僅與目前時段有關，還與其他時段有關。是以，需求彈性系數分為自彈性系數和交叉彈性系數，可表示為：

式中，

、

分别為

時段的自彈性系數及

時段對

時段的交叉彈性系數﹔

、

分别為原系統的負荷量和電價;

、

分别為系統負荷量和電價的改變量。在采用需求側響應後,各階段

負荷變化量

為，如下式：              

式中，下角标

、

、

代表伸縮系數歸屬的不同時間段。綜上。采用需求側響應之後各個時間的負荷需求量矩陣為

是以基于價格型需求響應成本為：

式中，

 為分時電價，

低谷電價為，

平段電價為，

高峰電價為.

3.2.3 冷能供應不足懲罰成本

為了避免對負荷冷能供應不足。可以施加懲罰項，這樣在所有階段都能盡可能滿足冷負荷的需求，對提供的冷能與所需的冷能的偏差施加二次懲罰 :

這裡

，未滿足冷負荷的價格($/(kWh

)，被模組化為與第 k 階段的建築物占用率成正比，與其熱品質成反比.

3.2.4 激勵性需求響應成本

可控負荷為一些用電時段具有一定靈活性的用電裝置,如洗衣機、可定時電鍋等，還有本文的制冰機。在排程時段内采用激勵型需求響應直接控制方式進行負荷轉移,并給予補償。激勵型需求響應成本為：

式中, 

($/kWh)為在激勵型需求響應計劃中，參與激勵型需求響應者獲得的價格,

是由需求響應計劃計算的第

階段的有功功率損耗。

淨階段成本是電能成本和冷能供應不足懲罰成本的總和:

需求響應成本是基于價格型需求響應成本和基于激勵型需求響應成本的總和: 

3.3 随機最優控制問題

給定初始狀态x0，特殊政策π={μ0，...，

}的總預期成本為:(其中u(K)=μk(x(K)))   

其中，對于所有 k ∈T，期望都超過了 w(k)。是以，随機最優控制問題為：

 其中限制适用于所有

.

4 算例及運作結果 

5 Matlab代碼實作 

%% 狀态和控制時間序列(圖 3. MPC 政策下的狀态和控制)

%====(1)繪制 x1, 冰罐的充電狀态=====
ks = 0 : 1 : N;
figure(iFigure); clf
subplot(5,1,1)
bar(ks,x(1,:),'r')
title('冰庫中的冷量', 'FontSize', titleSize)
ylabel('x_1(kWh_{cool})')
set(gca,'XLim',[0,N],'XTick',0:12:N,'XTickLabel',0:2:round(N/2))

%====(2)繪制x2, 為建築冷卻不足(kWhcool）即前階段所有未滿足的冷卻負荷的總和（滞後冷卻負荷、積壓的冷負荷。）
subplot(5,1,2)
bar(ks,x(2,:),'g') %繪制階梯圖
title('冷量供應不足', 'FontSize', titleSize)
ylabel('x_2(kWh_{cool})')
set(gca,'XLim',[0,N],'XTick',0:12:N,'XTickLabel',0:2:round(N/2))

%=====(3)繪制u1, 制冰機消耗的功率。======
subplot(5,1,3)
bar(ks(1:end-1),u(1,:),'b')
title('制冰機消耗的有功功率', 'FontSize', titleSize)
ylabel('u_1(kW)')
set(gca,'XLim',[0,N],'XTick',0:12:N,'XTickLabel',0:2:round(N/2))

%====(4)繪制u2, 冰融化提供的冷卻。=====
subplot(5,1,4)
bar(ks(1:end-1),u(2,:),'c')
title('融化冰所提供的冷量', 'FontSize', titleSize)
ylabel('u_2(kW_{cool})')
set(gca,'XLim',[0,N],'XTick',0:12:N,'XTickLabel',0:2:round(N/2))

%======繪制u3, 主制冰機消耗的功率。=====
subplot(5,1,5)
bar(ks(1:end-1),u(3,:),'m')
title('主制冰機消耗的有功功率', 'FontSize', titleSize)
ylabel('u_3(kW)')
set(gca,'XLim',[0,N],'XTick',0:12:N,'XTickLabel',0:2:round(N/2))
xlabel('時間(小時)')

%% 電力和成本時間序列；成本條形圖（圖 4. MPC 政策下的電力和成本）
%該子產品繪制了耗電量及其基線，以及階段成本和總成本條形圖。

%==(1)繪制總電力、其基線，以及如果所有冷負荷直接由主冷水機組（稱為pNaive）滿足，所消耗的有功功率。===

% pNaive = w(1,:)./params.kMain + w(2,:);
% figure(iFigure+1); clf
% subplot(2,1,1)
% hold on
% plot(ks(1:end-1), u(1,:) + u(3,:) + w(2,:), 'c.-')  %MPC
% plot(ks(1:end-1), params.pBase, 'b')  %OLOC方法,不考慮DR
% %plot(ks(1:end-1), pNaive, 'r')   %pNaive = (w(1,:)./TSDRparams.kMain + w(2,:));如果所有冷負荷直接由主冷水機組（稱為pNaive）滿足
% set(gca,'XLim',[0,N],'XTick',0:12:N,'XTickLabel',0:6:round(N/2), ...
%     'YLim',[0,50*ceil(max(pNaive)/50)])
% title('有功功率消耗', 'FontSize', titleSize)
% ylabel('有功功率（kW）')
% legend(legLabel,'随機開環最優控制','理想情況', ...   %,, '随機開環最優控制', '冷負荷全部由制冰機供冷',
%     'Location', 'SouthEast', 'Orientation', 'Horizontal')
% box on

%

%======(2)階段成本和總成本條形圖。==========

%  figure(iFigure+2); clf
figure(iFigure+1);
%subplot(2,1,2)
barsToPlot = [costs.energyCosts', costs.DRcosts', costs.disutilityCosts'];
bar(0:N-1, barsToPlot, 'stack')
title('成本', 'FontSize', titleSize)
xlabel('時間(小時)')
ylabel('成本')
legend('電能成本', '激勵型需求響應成本', '冷量供應不足懲罰成本', ...  
    'Location', 'SouthWest', 'Orientation', 'Horizontal')
set(gca, 'XLim', [0,N], 'XTick',0:12:N,'XTickLabel',0:6:round(N/2))