COMP2823 1 S1 2023R

Java Python COMP2823

Assignment 1

S1 2023

Problem 1. (20 points)

Given an array A consisting of n integers, we want to compute a matrix B where for any 0 ≤ i < j < n we have

B[i][j] = f([A[i], A[i + 1], ..., A[j − 1]])

Consider the following algorithm for computing B:

Algorithm 1 Range Function Computation

1: function RangeFunc(A)

2: B ← new n × n matrix

3: for i ← 0 to n − 1 do

4: for j ← i + 1 to n − 1 do

5: C ← make a copy of A[i : j]

6: B[i][j] ← f(C)

7: return B

Assume that f(C) runs in Θ(log |C|) time and that allocating the space for the matrix takes O(n 2 ) time.

Your task is to

a) Using O-notation, upperbound the running time of RangeFunc. Explain your answer with a detailed line by line analysis.

b) Using Ω-notation, lowerbound the running time of RangeFunc. Explain your answer.

Problem 2. (40 points)

We would like to design an augmented queue data structure. In addition to the usual enqueue and dequeue operations, you need to support the even-diff operation, which when run on a queue Q = ⟨q0, q1, q2, . . . , qn−1⟩ returns

Examples:

• even-diff([1, 3, 50, 48]) returns 4,

• even-diff([1, 3, 50, 48, 30]) returns 4,

• even-diff([3, 50, 48, 30]) returns 65.

We are to design an implementation of the methods enqueue, dequeue, and even-diff so that all operations run in O(1) time. You can assume that the data structure always starts from the empty queue.

Your data structure should take O(n) space, where n is the number of elements currently stored in the data structure.

Your task is to:

a) Design a data structure that supports the required operations in the required time and space.

b) Briefly argue the correctness of your data structure and operations.

c) Analyse the running time of your operations and space of your data structure.

Problem 3. (40 points)

We would like to implement an ADT for keeping track of a collection of queues.

We would like to implement the following operations:

• init: Initialize the system having a single empty queue

• enqueue(Q, e): Q is a position specifying a given queue in the system, and e is the element to be enqueued into Q

• dequeue(Q, e): Q is a position specifying a given queue in the system, from which we could like to dequeue an element

• split(Q): Q is a position specifying a given queue in the system whose elements needs to be redistributed evenly into two new queues.

• iterate(): iterate over the queues stored in the system

Examples:

• given the state [⟨1, 5, 4⟩,⟨3, 7, 5, 2⟩,⟨10⟩] if we split the second queue the state should become [⟨1, 5, 4⟩,⟨3, 7⟩,⟨5, 2⟩,⟨10⟩],

• dai 写COMP2823 Assignment 1 S1 2023R given [⟨1, 5, 4⟩,⟨3, 7, 5, 2⟩] if we dequeue from the first queue the state should become [⟨5, 4⟩,⟨3, 7, 5, 2⟩],

Design a data structures for the ADT such that enqueue, dequeue, and init run in O(1) time, and split runs in O(log m) amortised time, where m is the number of operations performed. Your data structure should take O(n) space, where n is the number of elements currently stored in the data structure.

Your task is to:

a) Design a data structure that supports the required operations in the required time and space.

b) Briefly argue the correctness of your data structure and operations.

c) Analyse the running time of your operations and space of your data structure.

Written Assignment Guidelines

• Assignments should be typed and submitted as pdf (no pdf containing text as images, no handwriting).

• Start by typing your student ID at the top of the first page of your submission. Do not type your name.

• Submit only your answers to the questions. Do not copy the questions.

• When asked to give a plain English description, describe your algorithm as you would to a friend over the phone, such that you completely and unambiguously describe your algorithm, including all the important (i.e., non-trivial) details. It often helps to give a very short (1-2 sentence) description of the overall idea, then to describe each step in detail. At the end you can also include pseudocode, but this is optional.

• In particular, when designing an algorithm or data structure, it might help you (and us) if you briefly describe your general idea, and after that you might want to develop and elaborate on details. If we don’t see/understand your general idea, we cannot give you marks for it.

• Be careful with giving multiple or alternative answers. If you give multiple answers, then we will give you marks only for "your worst answer", as this indicates how well you understood the question.

• Some of the questions are very easy (with the help of the slides or book). You can use the material presented in the lecture or book without proving it. You do not need to write more than necessary (see comment above).

• When giving answers to questions, always prove/explain/motivate your answers.

• When giving an algorithm as an answer, the algorithm does not have to be given as (pseudo-)code.

• If you do give (pseudo-)code, then you still have to explain your code and your ideas in plain English.

• Unless otherwise stated, we always ask about worst-case analysis, worst case running times, etc.

• As done in the lecture, and as it is typical for an algorithms course, we are interested in the most efficient algorithms and data structures.

• If you use further resources (books, scientific papers, the internet,...) to formulate your answers, then add references to your sources and explain it in your own words. Only citing a source doesn’t show your understanding and will thus get you very few (if any) marks. Copying from any source without reference is considered plagiarism         

请你检查以下代码,说明需要优化或修改的地方:%% 考虑零碳排放的综合能源系统优化调度方法(电网+热网) clc close %% 一、控制参数设置 N_Time = 24; % 调度周期(小时) S_base = 10; % 基准功率 MW V_base = 12.66; % 基准电压 kV Z_base = V_base^2/S_base; % 基准阻抗 %% ==================== 风电机组参数配置 ==================== Ind_gen = [2, 7, 19, 26]; % 风机并网的节点 P_WT_g1 = [6.88 7.08 7.20 7.16 6.96 6.52 6.44 5.98 5.72 5.54 5.36 5.12 ... 4.64 4.56 4.60 4.64 4.52 4.52 4.92 5.40 5.96 6.56 6.68 6.72]-4.2; % 风电机组 #1(MW) MW 3MW P_WT_g2 = [16.6 16.4 16.5 16.6 16.8 11.7 11.3 11.3 12.3 13.5 14.9 16.4 17.2 17.7 18 17.9 17.4 ... 16.3 16.1 16.2 16.6 16.8 16.9 16.8]/40; % 风电机组 #2(MW) 0.6 MW P_WT_g3 = P_WT_g2; P_WT_g4 = P_WT_g2; P_WT_g1 = P_WT_g1/S_base; % 标幺化处理 P_WT_g2 = P_WT_g2/S_base; P_WT_g3 = P_WT_g3/S_base; P_WT_g4 = P_WT_g4/S_base; P_WT_g = [P_WT_g1;P_WT_g2;P_WT_g3;P_WT_g4]; % 4台风机 %% ==================== 电网节点参数配置 ==================== % 电源母线参数 (MW、MVar) % No.(1) |Type(2)| Pd(3)| Qd(4) power_bus = [ 1 3 0 0; 2 1 0.100 0.060; 3 1 0.090 0.040; 4 1 0.120 0.080; 5 1 0.060 0.030; 6 1 0.060 0.020; 7 1 0.200 0.100; 8 1 0.200 0.100; 9 1 0.060 0.020; 10 1 0.060 0.020; 11 1 0.045 0.030; 12 1 0.060 0.035; 13 1 0.060 0.035; 14 1 0.120 0.080; 15 1 0.060 0.010; 16 1 0.060 0.020; 17 1 0.060 0.020; 18 1 0.090 0.040; 19 1 0.090 0.040; 20 1 0.090 0.040; 21 1 0.090 0.040; 22 1 0.090 0.040; 23 1 0.090 0.050; 24 1 0.420 0.200; 25 1 0.420 0.200; 26 1 0.060 0.025; 27 1 0.060 0.025; 28 1 0.060 0.020; 29 1 0.120 0.070; 30 1 0.200 0.100; 31 1 0.150 0.070; 32 1 0.210 0.100; 33 1 0.060 0.040; ]; N_bus_1 = size(power_bus,1); % 33节点网络 Pd_ratio = power_bus(:,3)/sum(power_bus(:,3)); % 有功负载比 Qd_ratio = power_bus(:,4)/sum(power_bus(:,4)); % 无功负载比 Pd_0 = [63 62 60 58 59 65 72 85 95 99 100 99 93 92 90 88 90 92 96 98 96 90 80 70]/15-1.5; % 基准有功功率 MW Qd_0 = [18 16 15 14 15.5 15 16 17 18 19 20 20.5 21 20.5 21 19.5 20 20 19.5 19.5 18.5 18.5 18 18]/10; % 基准无功功率 MVar Pd = Pd_ratio * Pd_0; % 节点电负荷计算 Qd = Qd_ratio * Qd_0; Pd = Pd/S_base; % 标幺化有功负载 Qd = Qd/S_base; % 标幺化无功负载 U2_min = 0.95^2; % 电压边界设置 U2_max = 1.05^2; %% ==================== 电网支路参数配置 ==================== % 电力线支路参数 % No.(1) |From bus(2 )|To bus(3) |r(4) |x(5) |P_line_max(6) |P_line_min(7) branch = [ 1 1 2 0.0922 0.0470 9.9 0; 2 2 3 0.4930 0.2512 9.9 0; 3 3 4 0.3661 0.1864 9.9 0; 4 4 5 0.3811 0.1941 9.9 0; 5 5 6 0.8190 0.7070 9.9 0; 6 6 7 0.1872 0.6188 9.9 0; 7 7 8 0.7115 0.2351 9.9 0; 8 8 9 1.0299 0.7400 9.9 0; 9 9 10 1.0440 0.7400 9.9 0; 10 10 11 0.1967 0.0651 9.9 0; 11 11 12 0.3744 0.1298 9.9 0; 12 12 13 1.4680 1.1549 9.9 0; 13 13 14 0.5416 0.7129 9.9 0; 14 14 15 0.5909 0.5260 9.9 0; 15 15 16 0.7462 0.5449 9.9 0; 16 16 17 1.2889 1.7210 9.9 0; 17 17 18 0.7320 0.5739 9.9 0; 18 2 19 0.1640 0.1565 9.9 0; 19 19 20 1.5042 1.3555 9.9 0; 20 20 21 0.4095 0.4784 9.9 0; 21 21 22 0.7089 0.9373 9.9 0; 22 3 23 0.4512 0.3084 9.9 0; 23 23 24 0.8980 0.7091 9.9 0; 24 24 25 0.8959 0.7071 9.9 0; 25 6 26 0.2031 0.1034 9.9 0; 26 26 27 0.2842 0.1447 9.9 0; 27 27 28 1.0589 0.9338 9.9 0; 28 28 29 0.8043 0.7006 9.9 0; 29 29 30 0.5074 0.2585 9.9 0; 30 30 31 0.9745 0.9629 9.9 0; 31 31 32 0.3105 0.3619 9.9 0; 32 32 33 0.3411 0.5302 9.9 0; ]; N_line = size(branch,1); B_line_i = branch(:,2); B_line_j = branch(:,3); r_pu = branch(:,4)/Z_base; % 电阻标幺值 x_pu = branch(:,5)/Z_base; % 电抗标幺值 P_max_pu = branch(:,6)/S_base; % 功率限值标幺化 P_min_pu = branch(:,7)/S_base; %% ==================== 无功补偿设备配置 ==================== % 电容无功补偿器参数 % Location(1) |Max(2) |Min(3) |Step(4) Com_Cap = [ 5 0.2 0 0.05; 10 0.2 0 0.05; 13 0.2 0 0.05; 17 0.2 0 0.05; 20 0.2 0 0.05; 23 0.2 0 0.05; 30 0.2 0 0.05; ]; v = 2; N_Com_Cap = size(Com_Cap,1); % 补偿器数量 Ind_Com_Cap = Com_Cap(:,1); % 补偿器位置 Com_Cap_Step = Com_Cap(:,4); % 电容步长调整 % 静止无功补偿器参数 SVG = [ % Location(1) |Max(2)| Min(3) 4 0.1 0; 9 0.1 0; 14 0.1 0; ]; Ind_SVG = SVG(:,1); % SVG数量 SVG(:, 2:3) = SVG(:, 2:3) / S_base; % 标幺化 %% ==================== 可调变压器配置 ==================== % 可调变压器参数 % Line No.(1) |K_max(2) |K_min(3) |K_Step(4)| OLTC = [ 1 1.05 0.95 0.01; 18 1.05 0.95 0.01; 22 1.05 0.95 0.01; 25 1.05 0.95 0.01; ]; t_OLTC = [0.95:0.01:1.05]; % 可用抽头值 T_OLTC = repmat(t_OLTC',1,N_Time); n_OLTC= length(t_OLTC); % 总损失价值 N_OLTC = size(OLTC,1); % 变压器数量 Ind_OLTC = OLTC(:,1); % 变压器所在线路索引 Ind_sub_line = zeros(N_line,2); % 电网子线路索引 for i = 1: N_line temp = find(B_line_i == B_line_j(i)); if ~isempty(temp) Ind_sub_line(i,1:length(temp)) = temp; end end Pg_min = zeros(N_bus_1,N_Time); Pg_max = zeros(N_bus_1,N_Time); %% ==================== 热网节点参数配置 ==================== % 热网节点参数(kW、par、℃) % No(1) |Hd(2)| Pr_SR_min(3)| tao_S_max(4)| tao_S_min(5)| tao_R_max(6)| tao_R_min(7)| mass flow(8) heat_node = [ 1 0 50000 120 90 80 60 0; 2 0 50000 120 90 80 60 0; 3 0 50000 120 90 80 60 0; 4 0 50000 120 90 80 60 0; 5 250 50000 120 90 80 60 2; 6 250 50000 120 90 80 60 2; 7 250 50000 120 90 80 60 2; 8 500 50000 120 90 80 60 4; ]; N_bus_2 = size(heat_node,1); H_ratio = heat_node(:,2)/sum(heat_node(:,2)); H_hd_0 = [1250*ones(1,4), 1150*ones(1,4), 1000*ones(1,4), 800*ones(1,4), 1150*ones(1,4), 1250*ones(1,4)]; % 热负荷曲线 kW H_Hd = H_ratio * H_hd_0; % 节点热负荷计算 Nd_Hd = find(heat_node(:,2)>0); m_Hd = heat_node(Nd_Hd,8); % 热负荷质量流量比 tao_NS_max = heat_node(:,4); % 温度边界设置 tao_NS_min = heat_node(:,5); tao_NR_max = heat_node(:,6); % 温度边界设置 tao_NR_min = heat_node(:,7); %% ==================== 热力发电参数配置 ==================== % 热力发电参数(kW、kg/s) % Location(1) |Hg_max(2) |Hg_min(3) |C_A(4) |C_B(5) |C_(6) |Mass flow(7) heat_gen = [ 2 2500 0 0.05 20 0 10 ]; N_gen = size(heat_gen,1); % 热力机组数量 Nd_HS = heat_gen(:,1); % 热力站的节点位置 mf_HS = heat_gen(:,7); % 热力机组质量流量 Hg_min = heat_gen(:, 3); % 最小出力 Hg_max = heat_gen(:, 2); % 最大出力 %% ==================== 热网管道参数配置 ==================== % 热网传输通道 % No.(1) |From node(2) |To node(3) |L(4) |u_p(5) |u_T(6) |ms_max(7) |ms_min(8) |mr_max(9) |mr_min(10) pipe = [ 1 1 2 1000 0.04 0.5*1e-3 1e6 0 10 10; 2 2 3 1000 0.04 0.5*1e-3 1e6 0 8 8; 3 3 4 1000 0.04 0.5*1e-3 1e6 0 6 6; 4 2 5 1000 0.04 0.5*1e-3 1e6 0 2 2; 5 3 6 1000 0.04 0.5*1e-3 1e6 0 2 2; 6 4 7 1000 0.04 0.5*1e-3 1e6 0 2 2; 7 4 8 1000 0.04 0.5*1e-3 1e6 0 4 4; ]; N_pipe = size(pipe,1); % 获取管道数量 pipe_i = pipe(:,2); % 提取节点i首端 pipe_j = pipe(:,3); % 提取节点j末端 L_pipe = pipe(:,4); % 管道长度 lamada_pipe = pipe(:,6); % 热传导系数 m_pipe_max = pipe(:, 7); % 最大流量 m_pipe_min = pipe(:, 8); % 最小流量 ms_pipe = pipe(:, 9); % 发送端流量 mr_pipe = pipe(:, 10); % 接收端流量 S_pipe_F = zeros(N_bus_2,2); % 发送端拓扑矩阵 S_pipe_T = zeros(N_bus_2,2); % 接收端拓扑矩阵 %% ==================== 压缩空气储能参数配置 ==================== % 集线器参数 CAES (kW) N_ct = 2; % 压缩机数量 N_et = 2; % 涡轮机数量 % 运行机组参数 yita_comp = [0.80, 0.75]; % 压缩机等熵效率 yita_turb = [0.86, 0.86]; % 涡轮机等熵效率 beta_comp = [11.6, 8.15]; % 压缩机压力比 pi_turb = [8.9, 8.9]; P_comp_min = zeros(N_ct, 1); P_comp_max = 500 * ones(N_ct, 1); P_turb_min = zeros(N_et, 1); P_turb_max = 1000 * ones(N_et, 1); % 物理常量参数 Vst = 2000; % m^3 k = 1.4; % 绝热指数 Rg = 0.297; % KJ/(kg.K) cp_a = 1.007; % KJ/(kg.K) cp_w = 4.2; % KJ/(kg.K) cp_s = 2.5; % KJ/(kg.K) % 环境温度参数 tao_K = 273.15; % 摄氏转开尔文偏移量 tao_am = 15 + tao_K; % 直接计算开尔文温度 tao_str = 40 + tao_K; % 介质初始温度 tao_salt_min = 60 + tao_K; tao_salt_max = 320 + tao_K; % 环境压力参数 pr_am = 0.101 * 1e3; % kPa pr_st_min = 8.4 * 1e3; % kPa pr_st_max = 9.0 * 1e3; % kPa % 质量流量参数 qm_comp_min = 0; qm_comp_max = 2.306 / 3.6; qm_turb_min = 0; qm_turb_max = 8.869 / 3.6; % 热功率上、下限参数 H_str_min = 0.2*1e3; % kW H_str_max = 3.0*1e3; % kW % ---------- 压力计算(向量化)---------- % 预压缩计算 k_ratio = (k - 1) / k; % 复用指数计算 % 压缩机压力计算 pr_comp_in_1 = pr_am * ones(1, N_Time); pr_comp_out_1 = beta_comp(1) * pr_comp_in_1; pr_comp_in_2 = pr_comp_out_1; pr_comp_out_2 = beta_comp(2) * pr_comp_in_2; % 涡轮机压力计算 pr_turb_in_1 = pr_st_min * ones(1, N_Time); pr_turb_out_1 = pr_turb_in_1 / pi_turb(1); pr_turb_in_2 = pr_turb_out_1; pr_turb_out_2 = pr_turb_in_2 / pi_turb(2); % ---------- 温度计算(复用计算因子)---------- % 压缩过程温度因子计算 y_comp = beta_comp .^ k_ratio; % 向量化计算 tao_comp_in_1 = tao_am * ones(1, N_Time); tao_comp_in_2 = (40 + tao_K) * ones(1, N_Time); tao_comp_out_1 = tao_comp_in_1 ./ yita_comp(1) .* (y_comp(1) - 1 + yita_comp(1)); tao_comp_out_2 = tao_comp_in_2 ./ yita_comp(2) .* (y_comp(2) - 1 + yita_comp(2)); % 热交换温度计算 tao_cold_s_in_1 = tao_comp_out_1; tao_cold_s_out_1 = 90 + tao_K; % 直接计算开尔文温度 tao_cold_w_in_1 = tao_cold_s_out_1; tao_cold_w_out_1 = 40 + tao_K; tao_cold_s_in_2 = tao_comp_out_2; tao_cold_s_out_2 = 90 + tao_K; tao_cold_w_in_2 = tao_cold_s_out_2; tao_cold_w_out_2 = 40 + tao_K; % 涡轮温度计算(复用计算因子) y_turb = pi_turb .^ (-k_ratio); % 向量化计算 tao_turb_in_1 = (280 + tao_K) * ones(1, N_Time); tao_turb_in_2 = (280 + tao_K) * ones(1, N_Time); tao_turb_out_1 = tao_turb_in_1 .* yita_turb(1) .* (y_turb(1) - 1 + 1./yita_turb(1)); tao_turb_out_2 = tao_turb_in_2 .* yita_turb(2) .* (y_turb(2) - 1 + 1./yita_turb(2)); % 热能传递计算 tao_heat_in_1 = tao_str; tao_heat_in_2 = tao_turb_out_1; tao_heat_out_1 = tao_turb_in_1; tao_heat_out_2 = tao_turb_in_2; % 系统标识 CAES_ind = 2; %% 二、运行变量设置 %% ==================== 电网运行变量 ==================== P_line_load = sdpvar(N_line,N_Time,'full'); % 每条线路的有功功率 Q_line_load = sdpvar(N_line,N_Time,'full'); % 每条线路的无功功率 P_sub = sdpvar(N_line,N_Time,'full'); Q_sub = sdpvar(N_line,N_Time,'full'); U2_node = sdpvar(N_bus_1,N_Time,'full'); % 母线电压幅值的平方 Pg_node = sdpvar(N_bus_1,N_Time,'full'); % 向每条母线注入发电机有功功率 Qc_node = sdpvar(N_bus_1,N_Time,'full'); % 无功补偿器SVG P_grid = sdpvar(1,N_Time); % 上网购电有功 Q_grid = sdpvar(1,N_Time); % 上网购电无功 %% ==================== 热网运行变量 ==================== tao_PS_src = sdpvar(N_pipe,N_Time,'full'); % 供水管“源”侧温度 tao_PS_load = sdpvar(N_pipe,N_Time,'full'); % 供水管“荷”侧温度 tao_PR_src = sdpvar(N_pipe,N_Time,'full'); % 回流管“源”侧温度 tao_PR_load = sdpvar(N_pipe,N_Time,'full'); % 回流管“荷”侧温度 tao_NS_node = sdpvar(N_bus_2,N_Time,'full'); % 供电网络节点温度 tao_NR_node = sdpvar(N_bus_2,N_Time,'full'); % 回流网络节点温度 Hg_HP = sdpvar(N_bus_2,N_Time,'full'); % 配备CAES的热泵的热功率 %% ==================== CAES系统运行变量 ==================== on_comp = binvar(1,N_Time); % 压缩机 comp 的开/关 on_turb = binvar(1,N_Time); % 涡轮机 turb 的开/关 P_comp_1 = sdpvar(1,N_Time); % 压缩机功率 comp1 P_comp_2 = sdpvar(1,N_Time); % 压缩机功率 comp2 P_turb_1 = sdpvar(1,N_Time); % 涡轮机功率 turb1 P_turb_2 = sdpvar(1,N_Time); % 涡轮机功率 turb2 P_caes_d = sdpvar(1,N_Time); P_caes_g = sdpvar(1,N_Time); %% ==================== 线性化辅助变量 ==================== pr_st = sdpvar(1,N_Time); % 储气罐压强 pr_st_0 = sdpvar(1,1); qm_comp = sdpvar(1,N_Time); qm_turb = sdpvar(1,N_Time); y1 = sdpvar(1,N_Time); % 线性化储气室压强约束 y2 = sdpvar(1,N_Time); h1 = sdpvar(1,N_Time); % 线性化储热系统SOC h2 = sdpvar(1,N_Time); HM = 1e7; % 大M法常数 H_coll_s1 = sdpvar(1,N_Time); H_coll_s2 = sdpvar(1,N_Time); H_cons1 = sdpvar(1,N_Time); H_cons2 = sdpvar(1,N_Time); H_coll_sum = sdpvar(1,N_Time); H_cons_sum = sdpvar(1,N_Time); H_str = sdpvar(1,N_Time); % 储热罐中存贮的热量 H_str_0 = sdpvar(1,1); Hg_CAES = sdpvar(1,N_Time); % 蓄热环节可供热负荷 % 无功补偿器分段线性化 for i = 1:N_Com_Cap % 电容器 xd_cap{i} = binvar(v+1,N_Time); delta_cap{i} = sdpvar(v+1,N_Time); end for i = 1:N_OLTC % 变压器 rd_tap{i} = binvar(n_OLTC,N_Time); h_tap{i} = sdpvar(n_OLTC,N_Time); end %% 三、约束条件设置 %% ==================== CAES系统约束 ==================== F_turb = []; % 每级功率定义 F_comp = []; % 每级功率定义 F_oper = []; % 运行约束 F_power = []; % 功率平衡约束 F_air_str = []; % 储气罐压强动态约束 F_cold = []; F_heat = []; F_heat_str = []; % 压缩机约束 F_comp = [F_comp, P_comp_1 == 1/yita_comp(1)*k/(k-1)*Rg*qm_comp.*tao_comp_in_1.*(y_comp1-1)]; F_comp = [F_comp, P_comp_2 == 1/yita_comp(2)*k/(k-1)*Rg*qm_comp.*tao_comp_in_2.*(y_comp2-1)]; F_comp = [F_comp, P_comp_min(1)*on_comp <= P_comp_1 <= P_comp_max(1)*on_comp ]; % 每级消耗的功率约束 F_comp = [F_comp, P_comp_min(2)*on_comp <= P_comp_2 <= P_comp_max(2)*on_comp ]; F_comp = [F_comp, P_caes_d == P_comp_1 + P_comp_2]; % 总功率定义 % 涡轮机约束 F_turb = [F_turb, P_turb_1 == yita_turb(1)*k/(k-1)*Rg*qm_turb.*tao_turb_in_1.*(1-y_turb1)]; F_turb = [F_turb, P_turb_2 == yita_turb(2)*k/(k-1)*Rg*qm_turb.*tao_turb_in_2.*(1-y_turb2)]; F_turb = [F_turb, P_turb_min(1)*on_turb <= P_turb_1 <= P_turb_max(1)*on_turb]; % 每级发出功率约束 F_turb = [F_turb, P_turb_min(2)*on_turb <= P_turb_2 <= P_turb_max(2)*on_turb]; F_turb = [F_turb, P_caes_g == P_turb_1 + P_turb_2]; % 机组运行约束 F_oper = [F_oper, 0 <= on_comp + on_turb <= 1]; % 充放电不能同时进行 F_oper = [F_oper, qm_comp_min*on_comp <= qm_comp <= qm_comp_max*on_comp]; % 质量流量非否约束 F_oper = [F_oper, qm_turb_min*on_turb <= qm_turb <= qm_turb_max*on_turb]; % 边界约束和质量流量约束 F_air_str = [F_air_str, pr_st(1) == pr_st_0 + 1/Vst * Rg * tao_str * 3600*(y1(1)- y2(1))]; F_air_str = [F_air_str, pr_st(1,2:N_Time) == pr_st(1,1:N_Time-1) + 1/Vst * Rg * tao_str * 3600*(y1(2:N_Time)- y2(2:N_Time))]; F_air_str = [F_air_str, pr_st(1,N_Time) == pr_st_0]; F_air_str = [F_air_str, qm_comp_min*on_comp <= y1 <= qm_comp_max*on_comp]; F_air_str = [F_air_str, qm_turb_min*on_turb <= y2 <= qm_turb_max*on_turb]; F_air_str = [F_air_str, qm_comp_min*(1-on_comp) <= qm_comp - y1 <= qm_comp_max*(1-on_comp)]; F_air_str = [F_air_str, qm_turb_min*(1-on_turb) <= qm_turb - y2 <= qm_turb_max*(1-on_turb)]; F_air_str = [F_air_str, pr_st_min <= pr_st <= pr_st_max]; F_air_str = [F_air_str, pr_st_min <= pr_st_0 <= pr_st_max]; % 热回收系统约束 - 压缩热回收 F_cold = [F_cold, H_coll_s1 == cp_a*qm_comp.*(tao_cold_s_in_1 - tao_cold_s_out_1)]; % 回热系统热量约束,冷却环节 F_cold = [F_cold, H_coll_s2 == cp_a*qm_comp.*(tao_cold_s_in_2 - tao_cold_s_out_2)]; % 第1,2级导热油回收的热量 F_cold = [F_cold, H_coll_sum == H_coll_s1 + H_coll_s2]; % 只计及导热油收集的热量,回收的总热量 % 热回收系统约束 - 膨胀热消耗 F_heat = [F_heat, H_cons1 == cp_a*qm_turb.*(tao_heat_out_1 - tao_heat_in_1)]; % 第1-2级消耗的热量 F_heat = [F_heat, H_cons2 == cp_a*qm_turb.*(tao_heat_out_2 - tao_heat_in_2)]; F_heat = [F_heat, H_cons_sum == H_cons1 + H_cons2]; % 加热环节消耗的总热量 % 线性化约束和边界约束 F_heat_str = [F_heat_str, H_str(1,1) == H_str_0 + h1(1,1) - h2(1,1) - Hg_CAES(1,1)]; % 回热系统SOC% 高温储热罐中存贮的热量 F_heat_str = [F_heat_str, H_str(1,2:N_Time) == H_str(1,1:N_Time-1) + h1(1,2:N_Time) - h2(1,2:N_Time) - Hg_CAES(1,2:N_Time)]; F_heat_str = [F_heat_str, -HM*on_comp <= h1 <= HM*on_comp]; F_heat_str = [F_heat_str, -HM*on_turb <= h2 <= HM*on_turb]; F_heat_str = [F_heat_str, -HM*(1-on_comp) <= H_coll_sum - h1 <= HM*(1-on_comp)]; F_heat_str = [F_heat_str, -HM*(1-on_turb) <= H_cons_sum - h2 <= HM*(1-on_turb)]; F_heat_str = [F_heat_str, H_str_min <= H_str <= H_str_max]; % Box 约束 F_heat_str = [F_heat_str, H_str_min <= H_str_0 <= H_str_max]; F_heat_str = [F_heat_str, H_str(1,N_Time) == H_str_0]; %% ==================== 电力网络约束 ==================== F_P = []; F_Q = []; F_U = []; % 节点功率约束 for t = 1:N_Time for i = 1:N_bus_1 % 设置未放置发电机的母线节点的注入(P,Q)功率及a,b,c为 0, 发电机与SVG节点约束 if isempty(find(Ind_gen == i)) F_P = [F_P, Pg_node(i,t) == 0]; else Pg_max(i,:) = P_WT_g(find(Ind_gen == i),:); end % 判断是否装有SVG if isempty(find(Ind_SVG ==i)) % 未安装SVG F_Q = [F_Q, Qc_node(i,t) == 0]; else % 安装SVG temp = find(Ind_SVG == i); % 索引 F_Q = [F_Q, SVG(temp,3) <= Qc_node(i,t) <= SVG(temp,2)]; end end end % 线路子线路潮流 for t = 1:N_Time for i = 1:N_line num_temp = size(find(Ind_sub_line(i,:) == 0),2); if num_temp == 1 F_P = [F_P, P_sub(i,t) == P_line_load(Ind_sub_line(i,1),t)]; F_Q = [F_Q, Q_sub(i,t) == Q_line_load(Ind_sub_line(i,1),t)]; elseif num_temp == 2 F_P = [F_P, P_sub(i,t) == 0]; F_Q = [F_Q, Q_sub(i,t) == 0]; else F_P = [F_P, P_sub(i,t) == P_line_load(Ind_sub_line(i,1),t) + P_line_load(Ind_sub_line(i,2),t)]; F_Q = [F_Q, Q_sub(i,t) == Q_line_load(Ind_sub_line(i,1),t) + Q_line_load(Ind_sub_line(i,2),t)]; end end end % 初始化计数器 OTLC_count = zeros(1,N_Time); ComCap_count = zeros(1,N_Time); M = 1000; % 基准节点约束 F_P = [F_P, P_line_load(1,:) == P_grid(1,:)]; F_Q = [F_Q, Q_line_load(1,:) == 0]; F_U = [F_U, U2_node(1,:) == 1.05^2]; Vs1 = 1.05^2 * ones(1,N_Time); yite = 0.8; % 母线功率平衡与设备约束 for t = 1:N_Time for i = 1:N_line if B_line_j(i) == CAES_ind CAES_ind; F_P = [F_P, P_line_load(i,t) + Pg_node(B_line_j(i),t) + P_caes_g(t)/1e3/S_base == P_sub(i,t) + Pd(B_line_j(i),t) + P_caes_d(t)/1e3/S_base]; else F_P = [F_P, P_line_load(i,t) + Pg_node(B_line_j(i),t) == P_sub(i,t) + Pd(B_line_j(i),t)]; end % 判断是否装有并联补偿电容 if ~isempty(find(Ind_Com_Cap == B_line_j(i))) % 安装补偿电容 ComCap_count(t) = ComCap_count(t) + 1; F_Q = [F_Q, Q_line_load(i,t) + 0.5*(U2_node(B_line_j(i),t)*Com_Cap(ComCap_count(t),3) + Com_Cap_Step(ComCap_count(t))*(2^0*delta_cap{ComCap_count(t)}(1,t) + ... 2^1*delta_cap{ComCap_count(t)}(2,t) + 2^2*delta_cap{ComCap_count(t)}(2,t)))+ Qc_node(B_line_j(i),t) == Q_sub(i,t) + Qd(B_line_j(i),t)]; for m = 1:v+1 F_U = [F_U, U2_node(B_line_j(i),t) - M*(1-xd_cap{ComCap_count(t)}(m,t)) <= delta_cap{ComCap_count(t)}(m,t) <= U2_node(B_line_j(i),t) + M*(1-xd_cap{ComCap_count(t)}(m,t))]; F_U = [F_U, -M*xd_cap{ComCap_count(t)}(m,t) <= delta_cap{ComCap_count(t)}(m,t) <= M*xd_cap{ComCap_count(t)}(m,t)]; end F_U = [F_U, 0 <= 2^0*xd_cap{ComCap_count(t)}(1,t)+ 2^1*xd_cap{ComCap_count(t)}(2,t) + 2^2*xd_cap{ComCap_count(t)}(3,t) <= (Com_Cap(ComCap_count(t),2) - ... Com_Cap(ComCap_count(t),3))/Com_Cap(ComCap_count(t),4)]; else % 未安装补偿电容 F_Q = [F_Q, Q_line_load(i,t) + Qc_node(B_line_j(i),t) == Q_sub(i,t) + Qd(B_line_j(i),t)]; end if ~isempty(find(Ind_OLTC == i)) % 含OTLC的支路 OTLC_count(t) = OTLC_count(t)+1; F_U = [F_U, sum(h_tap{OTLC_count(t)}(:,t)./T_OLTC(:,t).^2,1) == U2_node(B_line_i(i),t)-(r_pu(i)*P_line_load(i,t)+x_pu(i)*Q_line_load(i,t))/Vs1(1,t)]; for k = 1:n_OLTC F_U = [F_U, -M*(1-rd_tap{OTLC_count(t)}(k,t)) + U2_node(B_line_j(i),t) <= h_tap{OTLC_count(t)}(k,t) <= U2_node(B_line_j(i),t) + M*(1-rd_tap{OTLC_count(t)}(k,t))]; F_U = [F_U, -M*rd_tap{OTLC_count(t)}(k,t) <= h_tap{OTLC_count(t)}(k,t)<= M*rd_tap{OTLC_count(t)}(k,t)]; end F_U = [F_U,sum(rd_tap{OTLC_count(t)}(:,t),1) == 1]; else % 不含OTLC的支路 F_U = [F_U, U2_node(B_line_j(i),t)== U2_node(B_line_i(i),t)-(r_pu(i)*P_line_load(i,t)+x_pu(i)*Q_line_load(i,t))/Vs1(1,t)]; end % 线路潮流约束 F_P = [F_P, P_min_pu(i) <= P_line_load(i,t) <= P_max_pu(i)]; end for i = 1:N_bus_1 F_P = [F_P,Pg_min(i,t) <= Pg_node(i,t) <= Pg_max(i,t)]; F_U = [F_U, U2_min <= U2_node(i,t) <= U2_max]; end end F_P = [F_P, P_grid >= 0]; %% ==================== 热力站约束 ==================== F_H = []; Count_HS = 0; Count_Hd = 0; for j = 1:N_bus_2 if ~isempty(find(Nd_HS == j)) % 若该节点设置供热机组 Count_HS = Count_HS + 1; F_H = [F_H, Hg_HP(j,:) + Hg_CAES(1,:) == cp_w*mf_HS(Count_HS)*(tao_NS_node(j,:) - tao_NR_node(j,:))]; % 增加CAES产生的热能 F_H = [F_H, Hg_min <= Hg_HP(j,:) <= Hg_max]; elseif ~isempty(find(Nd_Hd == j)) % 若该节点设置供热负荷 Count_Hd = Count_Hd + 1; F_H = [F_H, Hg_HP(j,:) == zeros(1,N_Time)]; F_H = [F_H, cp_w*m_Hd(Count_Hd,:)*(tao_NS_node(j,:) - tao_NR_node(j,:)) == H_Hd(j,:)]; else % 其他节点 F_H = [F_H, Hg_HP(j,:) == zeros(1,N_Time)]; end F_H = [F_H, tao_NS_min(j) <= tao_NS_node(j,:) <= tao_NS_max(j)]; F_H = [F_H, tao_NR_min(j) <= tao_NR_node(j,:)<= tao_NR_max(j)]; end F_H = [F_H, 0 <= Hg_CAES <= 0.2*H_str]; %% ==================== 供热网络约束 ==================== F_PH = []; for i = 1:N_bus_2 temp1 = find(pipe_i == i); % 管道首节点为i temp2 = find(pipe_j == i); % 管道末节点为i S_pipe_F(i,1:length(temp1)) = temp1; % 首节点为i的管道集合 S_pipe_T(i,1:length(temp2)) = temp2; % 末节点为i的管道集合 end for i = 1:N_bus_2 % 供水管道温度节点混合 num_temp1 = size(find(S_pipe_T(i,:) == 0),2); if num_temp1 == 1 % 以节点i末端的管道数为1 b = S_pipe_T(i,1); % 管道编号 F_PH = [F_PH, ms_pipe(b)* tao_PS_load(b,:) == ms_pipe(b) * tao_NS_node(i,:)]; F_PH = [F_PH, tao_PR_src(b,:) == tao_NR_node(i,:)]; elseif num_temp1 == 0 % 以节点i末端的管道数为2 b = S_pipe_T(i,:); % 管道编号 F_PH = [F_PH, ms_pipe(b(1))* tao_PS_load(b(1),:) + ms_pipe(b(2))* tao_PS_load(b(2),:) == (ms_pipe(b(1))+ms_pipe(b(2))) * tao_NS_node(i,:)]; F_PH = [F_PH, tao_PR_src(b(1),:) == tao_NR_node(i,:)]; F_PH = [F_PH, tao_PR_src(b(2),:) == tao_NR_node(i,:)]; else % 以节点i末端的管道数为0 ,首节点 end % 回水管道温度节点混合 num_temp2 = size(find(S_pipe_F(i,:) == 0),2); if num_temp2 == 1 % 以节点i首端的管道数为1 b = S_pipe_F(i,1); % 管道编号 F_PH = [F_PH, mr_pipe(b)* tao_PR_load(b,:) == mr_pipe(b) * tao_NR_node(i,:)]; F_PH = [F_PH, tao_PS_src(b,:) == tao_NS_node(i,:)]; elseif num_temp2 == 0 % 以节点i首端的管道数为2 b = S_pipe_F(i,:); % 管道编号 F_PH = [F_PH, mr_pipe(b(1))* tao_PR_load(b(1),:) + mr_pipe(b(2))* tao_PR_load(b(2),:) == (mr_pipe(b(1))+mr_pipe(b(2))) * tao_NR_node(i,:)]; F_PH = [F_PH, tao_PS_src(b(1),:) == tao_NS_node(i,:)]; F_PH = [F_PH, tao_PS_src(b(2),:) == tao_NS_node(i,:)]; else % 以节点i末端的管道数为0 末节点 end end for i = 1:N_pipe % 温度变化方程 F_PH = [F_PH, tao_PS_load(i,:) == (tao_PS_src(i,:) - (tao_am-tao_K)).*exp(-lamada_pipe(i)*L_pipe(i)/(cp_w*ms_pipe(i))) + (tao_am-tao_K)]; F_PH = [F_PH, tao_PR_load(i,:) == (tao_PR_src(i,:) - (tao_am-tao_K)).*exp(-lamada_pipe(i)*L_pipe(i)/(cp_w*mr_pipe(i))) + (tao_am-tao_K)]; end
09-12
请检查以下代码,指出需要优化或者修改的地方:%% 考虑零碳排放的综合能源系统优化调度方法(热网+电网) clc close %% 设置参数 NT = 24; % 总调度周期 Nc = 2; % 压缩机数量 Ne = 2; % 涡轮机数量 S_base = 10; % 基准功率 MW V_base = 12.66; % 基准电压 kV Z_base = V_base^2/S_base; % 基准阻抗 I_base = S_base/(sqrt(3)*V_base); % 基准电流kA %% 电源总线 bus % No.(1) |Type(2)| Pd(3)| Qd(4) powerbus = [ 1 3 0 0; % MW、MVar 2 1 0.100 0.060; 3 1 0.090 0.040; 4 1 0.120 0.080; 5 1 0.060 0.030; 6 1 0.060 0.020; 7 1 0.200 0.100; 8 1 0.200 0.100; 9 1 0.060 0.020; 10 1 0.060 0.020; 11 1 0.045 0.030; 12 1 0.060 0.035; 13 1 0.060 0.035; 14 1 0.120 0.080; 15 1 0.060 0.010; 16 1 0.060 0.020; 17 1 0.060 0.020; 18 1 0.090 0.040; 19 1 0.090 0.040; 20 1 0.090 0.040; 21 1 0.090 0.040; 22 1 0.090 0.040; 23 1 0.090 0.050; 24 1 0.420 0.200; 25 1 0.420 0.200; 26 1 0.060 0.025; 27 1 0.060 0.025; 28 1 0.060 0.020; 29 1 0.120 0.070; 30 1 0.200 0.100; 31 1 0.150 0.070; 32 1 0.210 0.100; 33 1 0.060 0.040; ]; N_bus1 = size(powerbus,1); Pd_ratio = powerbus(:,3)/sum(powerbus(:,3)); % 有功负载比 Qd_ratio = powerbus(:,4)/sum(powerbus(:,4)); % 无功负载比 Pd0 = [63 62 60 58 59 65 72 85 95 99 100 99 93 92 90 88 90 92 96 98 96 90 80 70]/15-1.5; % 有功MW Qd0 = [18 16 15 14 15.5 15 16 17 18 19 20 20.5 21 20.5 21 19.5 20 20 19.5 19.5 18.5 18.5 18 18]/10; % 无功MVar Pd = Pd_ratio * Pd0; Qd = Qd_ratio * Qd0; Pd = Pd/S_base; Qd = Qd/S_base; U2_min = 0.95^2; U2_max = 1.05^2; %% 补偿器 % Location(1) |Max(2) |Min(3) |Step(4) ComCap = [ 5 0.2 0 0.05; 10 0.2 0 0.05; 13 0.2 0 0.05; 17 0.2 0 0.05; 20 0.2 0 0.05; 23 0.2 0 0.05; 30 0.2 0 0.05; ]; v = 2; N_ComCap = size(ComCap,1); % 补偿器数量 Ind_ComCap = ComCap(:,1); S = ComCap(:,4); %% 静止无功补偿器 SVG = [ % Mvar % Location(1) |Max(2)| Min(3) 4 0.1 0; 9 0.1 0; 14 0.1 0; ]; Ind_SVG = SVG(:,1); SVG(:,2:3) = SVG(:,2:3)/S_base; %% 热网节点 % No(1) |Hd(2)| Pr_SR_min(3)| tao_S_max(4)| tao_S_min(5)| tao_R_max(6)| tao_R_min(7)| mass flow(8) heatnode = [ %kW %par %℃ 1 0 50000 120 90 80 60 0; 2 0 50000 120 90 80 60 0; 3 0 50000 120 90 80 60 0; 4 0 50000 120 90 80 60 0; 5 250 50000 120 90 80 60 2; 6 250 50000 120 90 80 60 2; 7 250 50000 120 90 80 60 2; 8 500 50000 120 90 80 60 4; ]; N_bus2 = size(heatnode,1); H_ratio = heatnode(:,2)/sum(heatnode(:,2)); H_hd0 = [1250*ones(1,4), 1150*ones(1,4), 1000*ones(1,4), 800*ones(1,4), 1150*ones(1,4), 1250*ones(1,4)]; % kW H_Hd = H_ratio * H_hd0; Nd_Hd = find(heatnode(:,2)>0); m_Hd = heatnode(Nd_Hd,8); % 热负荷质量流量比 tao_NS_max = heatnode(:,4); tao_NS_min = heatnode(:,5); tao_NR_max = heatnode(:,6); tao_NR_min = heatnode(:,7); %% 电力发电 Ind_gen = [2 7 19 26]; Wg1 = [6.88 7.08 7.20 7.16 6.96 6.52 6.44 5.98 5.72 5.54 5.36 5.12 ... 4.64 4.56 4.60 4.64 4.52 4.52 4.92 5.40 5.96 6.56 6.68 6.72]-4.2; % Wind Gen #1(MW) MW 3MW Wg2 = [16.6 16.4 16.5 16.6 16.8 11.7 11.3 11.3 12.3 13.5 14.9 16.4 17.2 17.7 18 17.9 17.4 ... 16.3 16.1 16.2 16.6 16.8 16.9 16.8]/40; % Wind Gen #2(MW) 0.6 MW Wg3 = Wg2; Wg4 = Wg2; Wg1 = Wg1/S_base; Wg2 = Wg2/S_base; Wg3 = Wg3/S_base; Wg4 = Wg4/S_base; Wg = [Wg1;Wg2;Wg3;Wg4]; % 4台风机 %% 热力发电 % Location(1) |Hg_max(2) |Hg_min(3) |C_A(4) |C_B(5) |C_(6) |Mass flow(7) % heatgen = [ % kW % kg/s 2 2500 0 0.05 20 0 10 ]; N_gen = size(heatgen,1); Nd_HS = heatgen(:,1); % 热力站节点 m_HS = heatgen(:,7); Hg_min = heatgen(:,3); Hg_max = heatgen(:,2); %% CAES集线器 yita_comp = [0.80, 0.75]; yita_turb = [0.86, 0.86]; beta_comp = [11.6,8.15]; pi_turb = [8.9,8.9]; Pcomp_min = zeros(Nc,1); Pcomp_max = 500*ones(Nc,1); %kW Pturb_min = zeros(Ne,1); Pturb_max = 1000*ones(Ne,1); %kW Vst = 2000; % m^3 k = 1.4; % 绝热指数 Rg = 0.297; % KJ/(kg.K) cp_a = 1.007; % KJ/(kg.K),25℃ cp_w = 4.2; % KJ/(kg.K),25℃ cp_s = 2.5; % KJ/(kg.K),25℃ tao_am = 15; % 环境温度 tao_K = 273.15; tao_am = tao_am + tao_K; tao_str = 40; tao_str = tao_str + tao_K; tao_salt_min = 60 + tao_K; tao_salt_max = 320 + tao_K; pr_am = 0.101*1e3; % 环境压力 pr_st_min = 8.4*1e3; % Kpa pr_st_max = 9.0*1e3; % Kpa qm_comp_min = 0; qm_comp_max = 2.306/3.6; % kg/s 1MW CAES qm_turb_min = 0; qm_turb_max = 8.869/3.6; % kg/s 1MW CAES H_str_min = 0.2*1e3; %kW H_str_max = 3.0*1e3; %kW pr_comp_in1 = pr_am*ones(1,NT); % 固定压力 pr_comp_out1 = beta_comp(1)*pr_comp_in1; pr_comp_in2 = pr_comp_out1; pr_comp_out2 = beta_comp(2)*pr_comp_in2; y_comp1 = (beta_comp(1))^((k-1)/k); y_comp2 = (beta_comp(2))^((k-1)/k); pr_turb_in1 = pr_st_min*ones(1,NT); pr_turb_out1 = pr_turb_in1/pi_turb(1); pr_turb_in2 = pr_turb_out1; pr_turb_out2 = pr_turb_in2/pi_turb(2); y_turb1 = (pi_turb(1))^(-(k-1)/k); y_turb2 = (pi_turb(2))^(-(k-1)/k); tao_comp_in1 = tao_am*ones(1,NT); % 固定压力 tao_comp_in2 = (40 + tao_K)*ones(1,NT); tao_comp_out1 = tao_comp_in1/yita_comp(1).*(y_comp1-1+yita_comp(1)); tao_comp_out2 = tao_comp_in2/yita_comp(2).*(y_comp2-1+yita_comp(2)); tao_cold_s_in1 = tao_comp_out1; tao_cold_s_out1 = 90 + tao_K ; % 固定盐热交换器输出温度,90℃ tao_cold_w_in1 = tao_cold_s_out1; tao_cold_w_out1 = 40 + tao_K ; % 水热交换器输出温度,40℃ tao_cold_s_in2 = tao_comp_out2; tao_cold_s_out2 = 90 + tao_K; tao_cold_w_in2 = tao_cold_s_out2; tao_cold_w_out2 = 40 + tao_K ; tao_turb_in1 = (280 + tao_K)*ones(1,NT); tao_turb_in2 = (280 + tao_K)*ones(1,NT); tao_turb_out1 = tao_turb_in1*yita_turb(1).*(y_turb1-1+1/yita_turb(1)); tao_turb_out2 = tao_turb_in2*yita_turb(2).*(y_turb2-1+1/yita_turb(2)); tao_heat_in1 = tao_str; tao_heat_in2 = tao_turb_out1; tao_heat_out1 = tao_turb_in1; tao_heat_out2 = tao_turb_in2; CAES_ind = 2; %% 电力线支路branch % No.(1) |From bus(2 )|To bus(3) |r(4) |x(5) |P_line_max(6) |P_line_min(7) % branch = [ 1 1 2 0.0922 0.0470 9.9 0; 2 2 3 0.4930 0.2512 9.9 0; 3 3 4 0.3661 0.1864 9.9 0; 4 4 5 0.3811 0.1941 9.9 0; 5 5 6 0.8190 0.7070 9.9 0; 6 6 7 0.1872 0.6188 9.9 0; 7 7 8 0.7115 0.2351 9.9 0; 8 8 9 1.0299 0.7400 9.9 0; 9 9 10 1.0440 0.7400 9.9 0; 10 10 11 0.1967 0.0651 9.9 0; 11 11 12 0.3744 0.1298 9.9 0; 12 12 13 1.4680 1.1549 9.9 0; 13 13 14 0.5416 0.7129 9.9 0; 14 14 15 0.5909 0.5260 9.9 0; 15 15 16 0.7462 0.5449 9.9 0; 16 16 17 1.2889 1.7210 9.9 0; 17 17 18 0.7320 0.5739 9.9 0; 18 2 19 0.1640 0.1565 9.9 0; 19 19 20 1.5042 1.3555 9.9 0; 20 20 21 0.4095 0.4784 9.9 0; 21 21 22 0.7089 0.9373 9.9 0; 22 3 23 0.4512 0.3084 9.9 0; 23 23 24 0.8980 0.7091 9.9 0; 24 24 25 0.8959 0.7071 9.9 0; 25 6 26 0.2031 0.1034 9.9 0; 26 26 27 0.2842 0.1447 9.9 0; 27 27 28 1.0589 0.9338 9.9 0; 28 28 29 0.8043 0.7006 9.9 0; 29 29 30 0.5074 0.2585 9.9 0; 30 30 31 0.9745 0.9629 9.9 0; 31 31 32 0.3105 0.3619 9.9 0; 32 32 33 0.3411 0.5302 9.9 0; ]; N_line = size(branch,1); line_i = branch(:,2); line_j = branch(:,3); r = branch(:,4)/Z_base; x = branch(:,5)/Z_base; Pmax = branch(:,6)/S_base; Pmin = branch(:,7)/S_base; %% 可调变压器 % Line No.(1) |K_max(2) |K_min(3) |K_Step(4)|% OLTC = [ 1 1.05 0.95 0.01; 18 1.05 0.95 0.01; 22 1.05 0.95 0.01; 25 1.05 0.95 0.01; ]; t_OLTC = 0.95:0.01:1.05; % 可用抽头值 T_OLTC = repmat(t_OLTC',1,NT); n_OLTC= length(t_OLTC); % 总损失价值 N_OLTC = size(OLTC,1);% 变压器数量 Ind_OLTC = OLTC(:,1); Ind_subline = zeros(N_line,2); % 电网子线路索引 for i = 1: N_line temp = find(line_i == line_j(i)); if ~isempty(temp) Ind_subline(i,1:length(temp)) = temp; end end Pg_min = zeros(N_bus1,NT); Pg_max = zeros(N_bus1,NT); %% 传输通道 % No.(1) |From node(2) |To node(3) |L(4) |u_p(5) |u_T(6) |ms_max(7) |ms_min(8) |mr_max(9) |mr_min(10) pipe = [ 1 1 2 1000 0.04 0.5*1e-3 1e6 0 10 10; 2 2 3 1000 0.04 0.5*1e-3 1e6 0 8 8; 3 3 4 1000 0.04 0.5*1e-3 1e6 0 6 6; 4 2 5 1000 0.04 0.5*1e-3 1e6 0 2 2; 5 3 6 1000 0.04 0.5*1e-3 1e6 0 2 2; 6 4 7 1000 0.04 0.5*1e-3 1e6 0 2 2; 7 4 8 1000 0.04 0.5*1e-3 1e6 0 4 4; ]; N_pipe = size(pipe,1); pipe_i = pipe(:,2); pipe_j = pipe(:,3); L_pipe = pipe(:,4); lamada_pipe = pipe(:,6); m_pipe_max = pipe(:,7); m_pipe_min = pipe(:,8); ms_pipe = pipe(:,9); mr_pipe = pipe(:,10); S_pipe_F = zeros(N_bus2,2); S_pipe_T = zeros(N_bus2,2); %% 运行变量 % % PDN P = sdpvar(N_line,NT,'full'); % 每条线路的有功功率 Q = sdpvar(N_line,NT,'full'); % 每条线路的无功功率 Psub = sdpvar(N_line,NT,'full'); Qsub = sdpvar(N_line,NT,'full'); I2 = sdpvar(N_line,NT,'full'); % 线路电流幅值的平方 U2 = sdpvar(N_bus1,NT,'full'); % 母线电压幅值的平方 Pg = sdpvar(N_bus1,NT,'full'); % 向每条母线注入发电机有功功率 Qc = sdpvar(N_bus1,NT,'full'); % 无功补偿器SVG Pgrid = sdpvar(1,NT); % 上网购电 Qgrid = sdpvar(1,NT); % % DHN tao_PS_F = sdpvar(N_pipe,NT,'full'); % 供水管“源”侧温度 tao_PS_T = sdpvar(N_pipe,NT,'full'); %供水管“荷”侧温度 tao_PR_F = sdpvar(N_pipe,NT,'full'); %回流管“源”侧温度 tao_PR_T = sdpvar(N_pipe,NT,'full'); %回流管“荷”侧温度 tao_NS = sdpvar(N_bus2,NT,'full'); % 供电网络节点温度 tao_NR = sdpvar(N_bus2,NT,'full'); % 回流网络节点温度 Hg_HP = sdpvar(N_bus2,NT,'full'); % 配备CAES的热泵的热功率 % % CAES on_comp = binvar(1,NT); % comp 的开/关 on_turb = binvar(1,NT); % turb 的开/关 Pcomp1 = sdpvar(1,NT); % 电源缺点 comp1 Pcomp2 = sdpvar(1,NT); Pturb1= sdpvar(1,NT); Pturb2= sdpvar(1,NT); Pcaes_d = sdpvar(1,NT); Pcaes_g = sdpvar(1,NT); pr_st = sdpvar(1,NT); % 储气罐压强 pr_st0 = sdpvar(1,1); qm_comp = sdpvar(1,NT); qm_turb = sdpvar(1,NT); y1 = sdpvar(1,NT); % 线性化储气室压强约束 y2 = sdpvar(1,NT); h1 = sdpvar(1,NT); % 线性化储热系统SOC h2 = sdpvar(1,NT); HM = 1e7; % big M H_coll_s1 = sdpvar(1,NT); H_coll_s2 = sdpvar(1,NT); H_cons1 = sdpvar(1,NT); H_cons2 = sdpvar(1,NT); H_coll_sum = sdpvar(1,NT); H_cons_sum = sdpvar(1,NT); H_str = sdpvar(1,NT); % 储热罐中存贮的热量 H_str0 = sdpvar(1,1); Hg_CAES = sdpvar(1,NT); % 蓄热环节可供热负荷 for i = 1:N_ComCap %% 线性化变量 xd{i} = binvar(v+1,NT); delta{i} = sdpvar(v+1,NT); end for i = 1:N_OLTC rd{i} = binvar(n_OLTC,NT); h{i} = sdpvar(n_OLTC,NT); end %% 约束条件 % % CAES F_turb = []; % 每级功率定义 F_comp = []; % 每级功率定义 F_oper = []; % 运行约束 F_power = []; % 功率平衡约束 F_airstr = []; % 储气罐压强动态约束 F_cold = []; F_heat = []; F_heatstr = []; F_comp = [F_comp, Pcomp1 == 1/yita_comp(1)*k/(k-1)*Rg*qm_comp.*tao_comp_in1.*(y_comp1-1)]; F_comp = [F_comp, Pcomp2 == 1/yita_comp(2)*k/(k-1)*Rg*qm_comp.*tao_comp_in2.*(y_comp2-1)]; F_comp = [F_comp, Pcomp_min(1)*on_comp <= Pcomp1 <= Pcomp_max(1)*on_comp ]; % 每级消耗的功率约束 F_comp = [F_comp, Pcomp_min(2)*on_comp <= Pcomp2 <= Pcomp_max(2)*on_comp ]; F_comp = [F_comp, Pcaes_d == Pcomp1 + Pcomp2]; % 总功率定义 F_turb = [F_turb, Pturb1 == yita_turb(1)*k/(k-1)*Rg*qm_turb.*tao_turb_in1.*(1-y_turb1)]; F_turb = [F_turb, Pturb2 == yita_turb(2)*k/(k-1)*Rg*qm_turb.*tao_turb_in2.*(1-y_turb2)]; F_turb = [F_turb, Pturb_min(1)*on_turb <= Pturb1 <= Pturb_max(1)*on_turb]; % 每级发出功率约束 F_turb = [F_turb, Pturb_min(2)*on_turb <= Pturb2 <= Pturb_max(2)*on_turb]; F_turb = [F_turb, Pcaes_g == Pturb1 + Pturb2]; F_oper = [F_oper, 0 <= on_comp + on_turb <= 1]; %充放电不能同时进行 F_oper = [F_oper, qm_comp_min*on_comp <= qm_comp <= qm_comp_max*on_comp]; %质量流量非否约束 F_oper = [F_oper, qm_turb_min*on_turb <= qm_turb <= qm_turb_max*on_turb]; F_airstr = [F_airstr, pr_st(1) == pr_st0 + 1/Vst * Rg * tao_str * 3600*(y1(1)- y2(1))]; F_airstr = [F_airstr, pr_st(1,2:NT) == pr_st(1,1:NT-1) + 1/Vst * Rg * tao_str * 3600*(y1(2:NT)- y2(2:NT))]; F_airstr = [F_airstr, pr_st(1,NT) == pr_st0]; F_airstr = [F_airstr, qm_comp_min*on_comp <= y1 <= qm_comp_max*on_comp]; F_airstr = [F_airstr, qm_turb_min*on_turb <= y2 <= qm_turb_max*on_turb]; F_airstr = [F_airstr, qm_comp_min*(1-on_comp) <= qm_comp - y1 <= qm_comp_max*(1-on_comp)]; F_airstr = [F_airstr, qm_turb_min*(1-on_turb) <= qm_turb - y2 <= qm_turb_max*(1-on_turb)]; F_airstr = [F_airstr, pr_st_min <= pr_st <= pr_st_max]; F_airstr = [F_airstr, pr_st_min <= pr_st0 <= pr_st_max]; F_cold = [F_cold, H_coll_s1 == cp_a*qm_comp.*(tao_cold_s_in1 - tao_cold_s_out1)]; % 回热系统热量约束,冷却环节 F_cold = [F_cold, H_coll_s2 == cp_a*qm_comp.*(tao_cold_s_in2 - tao_cold_s_out2)]; % 第1,2级导热油回收的热量 F_cold = [F_cold, H_coll_sum == H_coll_s1 + H_coll_s2]; % 只计及导热油收集的热量,回收的总热量 F_heat = [F_heat, H_cons1 == cp_a*qm_turb.*(tao_heat_out1 - tao_heat_in1)]; % 第1-2级消耗的热量 F_heat = [F_heat, H_cons2 == cp_a*qm_turb.*(tao_heat_out2 - tao_heat_in2)]; F_heat = [F_heat, H_cons_sum == H_cons1 + H_cons2]; % 加热环节消耗的总热量 F_heatstr = [F_heatstr, H_str(1,1) == H_str0 + h1(1,1) - h2(1,1) - Hg_CAES(1,1)]; % 回热系统SOC% 高温储热罐中存贮的热量 F_heatstr = [F_heatstr, H_str(1,2:NT) == H_str(1,1:NT-1) + h1(1,2:NT) - h2(1,2:NT) - Hg_CAES(1,2:NT)]; F_heatstr = [F_heatstr, -HM*on_comp <= h1 <= HM*on_comp]; F_heatstr = [F_heatstr, -HM*on_turb <= h2 <= HM*on_turb]; F_heatstr = [F_heatstr, -HM*(1-on_comp) <= H_coll_sum - h1 <= HM*(1-on_comp)]; F_heatstr = [F_heatstr, -HM*(1-on_turb) <= H_cons_sum - h2 <= HM*(1-on_turb)]; F_heatstr = [F_heatstr, H_str_min <= H_str <= H_str_max]; % Box 约束 F_heatstr = [F_heatstr, H_str_min <= H_str0 <= H_str_max]; F_heatstr = [F_heatstr, H_str(1,NT) == H_str0]; %% PDN 约束 F_P = []; F_Q = []; F_U = []; for t = 1:NT for i = 1:N_bus1 % 设置未放置发电机的母线节点的注入(P,Q)功率及a,b,c为 0 if isempty(find(Ind_gen == i)) F_P = [F_P, Pg(i,t) == 0]; else Pg_max(i,:) = Wg(find(Ind_gen == i),:); end % 判断是否装有SVG if isempty(find(Ind_SVG ==i)) % 未安装SVG F_Q = [F_Q, Qc(i,t) == 0]; else % 安装SVG temp = find(Ind_SVG == i); % 索引 F_Q = [F_Q, SVG(temp,3) <= Qc(i,t) <= SVG(temp,2)]; end end end % 线路子线路潮流 for t = 1:NT for i = 1:N_line num_temp = size(find(Ind_subline(i,:) == 0),2); if num_temp == 1 F_P = [F_P, Psub(i,t) == P(Ind_subline(i,1),t)]; F_Q = [F_Q, Qsub(i,t) == Q(Ind_subline(i,1),t)]; elseif num_temp == 2 F_P = [F_P, Psub(i,t) == 0]; F_Q = [F_Q, Qsub(i,t) == 0]; else F_P = [F_P, Psub(i,t) == P(Ind_subline(i,1),t) + P(Ind_subline(i,2),t)]; F_Q = [F_Q, Qsub(i,t) == Q(Ind_subline(i,1),t) + Q(Ind_subline(i,2),t)]; end end end OTLC_count = zeros(1,NT); ComCap_count = zeros(1,NT); M = 1000; F_P = [F_P, P(1,:) == Pgrid(1,:)]; F_Q = [F_Q, Q(1,:) == 0]; F_U = [F_U, U2(1,:) == 1.05^2]; Vs1 = 1.05^2 * ones(1,NT); yite = 0.8; for t = 1:NT for i = 1:N_line if line_j(i) == CAES_ind CAES_ind; F_P = [F_P, P(i,t) + Pg(line_j(i),t) + Pcaes_g(t)/1e3/S_base == Psub(i,t) + Pd(line_j(i),t) + Pcaes_d(t)/1e3/S_base]; else F_P = [F_P, P(i,t) + Pg(line_j(i),t) == Psub(i,t) + Pd(line_j(i),t)]; end % 判断是否装有并联补偿电容 if ~isempty(find(Ind_ComCap == line_j(i))) % 安装补偿电容 ComCap_count(t) = ComCap_count(t) + 1; F_Q = [F_Q, Q(i,t) + 0.5*(U2(line_j(i),t)*ComCap(ComCap_count(t),3) + S(ComCap_count(t))*(2^0*delta{ComCap_count(t)}(1,t) + ... 2^1*delta{ComCap_count(t)}(2,t) + 2^2*delta{ComCap_count(t)}(2,t)))+ Qc(line_j(i),t) == Qsub(i,t) + Qd(line_j(i),t)]; for m = 1:v+1 F_U = [F_U, U2(line_j(i),t) - M*(1-xd{ComCap_count(t)}(m,t)) <= delta{ComCap_count(t)}(m,t) <= U2(line_j(i),t) + M*(1-xd{ComCap_count(t)}(m,t))]; F_U = [F_U, -M*xd{ComCap_count(t)}(m,t) <= delta{ComCap_count(t)}(m,t) <= M*xd{ComCap_count(t)}(m,t)]; end F_U = [F_U, 0 <= 2^0*xd{ComCap_count(t)}(1,t)+ 2^1*xd{ComCap_count(t)}(2,t) + 2^2*xd{ComCap_count(t)}(3,t) <= (ComCap(ComCap_count(t),2) - ... ComCap(ComCap_count(t),3))/ComCap(ComCap_count(t),4)]; else % 未安装补偿电容 F_Q = [F_Q, Q(i,t) + Qc(line_j(i),t) == Qsub(i,t) + Qd(line_j(i),t)]; end if ~isempty(find(Ind_OLTC == i)) % 含OTLC的支路 OTLC_count(t) = OTLC_count(t)+1; F_U = [F_U, sum(h{OTLC_count(t)}(:,t)./T_OLTC(:,t).^2,1) == U2(line_i(i),t)-(r(i)*P(i,t)+x(i)*Q(i,t))/Vs1(1,t)]; for k = 1:n_OLTC F_U = [F_U, -M*(1-rd{OTLC_count(t)}(k,t)) + U2(line_j(i),t) <= h{OTLC_count(t)}(k,t) <= U2(line_j(i),t) + M*(1-rd{OTLC_count(t)}(k,t))]; F_U = [F_U, -M*rd{OTLC_count(t)}(k,t) <= h{OTLC_count(t)}(k,t)<= M*rd{OTLC_count(t)}(k,t)]; end F_U = [F_U,sum(rd{OTLC_count(t)}(:,t),1) == 1]; else % 不含OTLC的支路 F_U = [F_U, U2(line_j(i),t)== U2(line_i(i),t)-(r(i)*P(i,t)+x(i)*Q(i,t))/Vs1(1,t)]; end % 线路潮流约束 F_P = [F_P, Pmin(i) <= P(i,t) <= Pmax(i)]; end for i = 1:N_bus1 F_P = [F_P,Pg_min(i,t) <= Pg(i,t) <= Pg_max(i,t)]; F_U = [F_U, U2_min <= U2(i,t) <= U2_max]; end end F_P = [F_P, Pgrid >= 0]; % 热力站 F_H = []; Count_HS = 0; Count_Hd = 0; for j = 1:N_bus2 if ~isempty(find(Nd_HS == j)) % 若该节点设置供热机组 Count_HS = Count_HS + 1; F_H = [F_H, Hg_HP(j,:) + Hg_CAES(1,:) == cp_w*m_HS(Count_HS)*(tao_NS(j,:) - tao_NR(j,:))]; % 增加CAES产生的热能 F_H = [F_H, Hg_min <= Hg_HP(j,:) <= Hg_max]; elseif ~isempty(find(Nd_Hd == j)) % 若该节点设置供热负荷 Count_Hd = Count_Hd + 1; F_H = [F_H, Hg_HP(j,:) == zeros(1,NT)]; F_H = [F_H, cp_w*m_Hd(Count_Hd,:)*(tao_NS(j,:) - tao_NR(j,:)) == H_Hd(j,:)]; else % 其他节点 F_H = [F_H, Hg_HP(j,:) == zeros(1,NT)]; end F_H = [F_H, tao_NS_min(j) <= tao_NS(j,:) <= tao_NS_max(j)]; F_H = [F_H, tao_NR_min(j) <= tao_NR(j,:)<= tao_NR_max(j)]; end F_H = [F_H, 0 <= Hg_CAES <= 0.2*H_str]; % 供热网络 F_PH = []; for i = 1:N_bus2 temp1 = find(pipe_i == i); % 管道首节点为i temp2 = find(pipe_j == i); % 管道末节点为i S_pipe_F(i,1:length(temp1)) = temp1; % 首节点为i的管道集合 S_pipe_T(i,1:length(temp2)) = temp2; % 末节点为i的管道集合 end for i = 1:N_bus2 % 供水管道温度节点混合 num_temp1 = size(find(S_pipe_T(i,:) == 0),2); if num_temp1 == 1 % 以节点i末端的管道数为1 b = S_pipe_T(i,1); % 管道编号 F_PH = [F_PH, ms_pipe(b)* tao_PS_T(b,:) == ms_pipe(b) * tao_NS(i,:)]; F_PH = [F_PH, tao_PR_F(b,:) == tao_NR(i,:)]; elseif num_temp1 == 0 % 以节点i末端的管道数为2 b = S_pipe_T(i,:); % 管道编号 F_PH = [F_PH, ms_pipe(b(1))* tao_PS_T(b(1),:) + ms_pipe(b(2))* tao_PS_T(b(2),:) == (ms_pipe(b(1))+ms_pipe(b(2))) * tao_NS(i,:)]; F_PH = [F_PH, tao_PR_F(b(1),:) == tao_NR(i,:)]; F_PH = [F_PH, tao_PR_F(b(2),:) == tao_NR(i,:)]; else % 以节点i末端的管道数为0 ,首节点 end % 回水管道温度节点混合 num_temp2 = size(find(S_pipe_F(i,:) == 0),2); if num_temp2 == 1 % 以节点i首端的管道数为1 b = S_pipe_F(i,1); % 管道编号 F_PH = [F_PH, mr_pipe(b)* tao_PR_T(b,:) == mr_pipe(b) * tao_NR(i,:)]; F_PH = [F_PH, tao_PS_F(b,:) == tao_NS(i,:)]; elseif num_temp2 == 0 % 以节点i首端的管道数为2 b = S_pipe_F(i,:); % 管道编号 F_PH = [F_PH, mr_pipe(b(1))* tao_PR_T(b(1),:) + mr_pipe(b(2))* tao_PR_T(b(2),:) == (mr_pipe(b(1))+mr_pipe(b(2))) * tao_NR(i,:)]; F_PH = [F_PH, tao_PS_F(b(1),:) == tao_NS(i,:)]; F_PH = [F_PH, tao_PS_F(b(2),:) == tao_NS(i,:)]; else % 以节点i末端的管道数为0 末节点 end end for i = 1:N_pipe % 温度变化方程 F_PH = [F_PH, tao_PS_T(i,:) == (tao_PS_F(i,:) - (tao_am-tao_K)).*exp(-lamada_pipe(i)*L_pipe(i)/(cp_w*ms_pipe(i))) + (tao_am-tao_K)]; F_PH = [F_PH, tao_PR_T(i,:) == (tao_PR_F(i,:) - (tao_am-tao_K)).*exp(-lamada_pipe(i)*L_pipe(i)/(cp_w*mr_pipe(i))) + (tao_am-tao_K)]; end %% 目标函数 Obj = 0; % 发电成本最小 C_grid = [0.05*ones(1,8) 0.10*ones(1,6) 0.08*ones(1,8) 0.05*ones(1,2)]; C_grid = C_grid * 1e4; Obj1 = C_grid*Pgrid'; % PDN 购电成本 Obj2 = C_grid*[(Hg_HP(CAES_ind,:))/1e3/S_base/yite]'; Obj = Obj1 + Obj2;
09-12
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值