total_load函数

rspower/examples/test_total_load.txt

+#include ../data/case14.txt
+#include ../lib/idx_gen.txt
+#include ../lib/idx_bus.txt
+#include ../lib/idx_brch.txt
+#include ../lib/total_load.txt
+#include ../lib/make_sdzip.txt
+
+//Pd = total_load(bus, gen, 1, [], [], []);
+
+on_temp = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);
+index_gbus_type = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on_temp) - 1;
+on = find(get_multi(slice(bus, [0], [BUS_TYPE-1, BUS_TYPE]), index_gbus_type) != PQ)';
+gbus = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on)' - 1;
+Pd = total_load(select(bus, gbus), [], 1, [], [], []);
+
+assert_relative_eq(Pd, [
+    [c(0,0)],
+    [c(21.7,0)],
+    [c(94.2,0)],
+    [c(11.2,0)],
+    [c(0,0)]
+]);
+
+return Pd;
\ No newline at end of file

rspower/lib/total_load.txt

@@ -73,8 +73,112 @@
 //       Pd = total_load(mpc, load_zone, opt)
 //
 //   See also SCALE_LOAD.
+//  调用本函数时，需要调用如下文件
+//  #include ../lib/idx_gen.txt
+//  #include ../lib/idx_bus.txt
+//  #include ../lib/idx_brch.txt
+//  #include ../lib/make_sdzip.txt
 
-fn  total_load(bus, gen, load_zone, opt, mpopt){
+fn  total_load(bus, gen, load_zone, opt_type, opt_nominal, mpopt){
+    //matlab中load_zone、opt为字符串输入，本程序中规定直接用数字来指代相应的选项类型
+    //load_zone为分区模式，包括、整个系统为一个分区(all，默认值，用0表示)、每个母线单独分区(bus，用1表示)、使用母线区域定义(area，用2表示)
+    //opt_type为负荷类型，包括固定+可调度负荷(BOTH，默认值，用0表示)、仅计算固定负荷(FIXED，用1表示)、仅计算可调度负荷(DISPATCHABLE，用2表示)
+    //opt_nominal为实际值(默认值，用0表示)，额定值(用1表示)
+
+    nb = size(bus, 0);         // number of buses
+
+    // default options
+    if is_empty(opt_type){
+        if is_empty(gen){
+            opt_type = 1;
+        }else{
+            opt_type = 0;
+        }
+    }
+    if is_empty(opt_nominal){
+        opt_nominal = 0;
+    }
+
+    want_fixed  = (opt_type == 0 || opt_type == 1); //是否计算固定负荷
+    want_disp   = (opt_type == 0 || opt_type == 2); //是否计算可调度负荷
+
+    // initialize load_zone
+    if is_empty(load_zone){
+        LZone = slice(bus, [0], [BUS_AREA-1, BUS_AREA]);
+    }else{
+        if load_zone == 0{          // make a single zone of all buses
+            LZone = ones(nb, 1);
+        }else if load_zone == 1{     // each bus is its own zone
+            LZone = range(1, nb+1)';
+        }else if load_zone == 2{     // use areas defined in bus data as zones
+            LZone = slice(bus, [0], [BUS_AREA-1, BUS_AREA]);
+        }
+    }
+    
+
+    nz = max(LZone);    // number of load zones
+
+    // fixed load at each bus, & initialize dispatchable
+    if want_fixed == 1{
+        Sd = make_sdzip(1, bus);
+        Vm = slice(bus, [0], [VM-1, VM]);
+        Sbusd = slice(Sd, [0], [2, 3]) + slice(Sd, [0], [1, 2]) .* Vm + slice(Sd, [0], [0, 1]) .* Vm.^2;
+        Pdf = real(Sbusd);      // real power
+        Qdf = imag(Sbusd);      // reactive power
+    }else{
+        Pdf = zeros(nb, 1);     // real power
+        Qdf = zeros(nb, 1);     // reactive power
+    }
+
+    // dispatchable load at each bus 
+    if want_disp ==1{           // need dispatchable
+        ng = size(gen, 0);
+        is_ld = (slice(gen, [0], [PMIN-1, PMIN]) < 0) && (slice(gen, [0], [PMAX-1, PMAX]) == 0) && (slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);
+        ld = find(is_ld);
+
+        // create map of external bus numbers to bus indices(创建母线编号映射)
+        i2e = slice(bus, [0], [BUS_I-1, BUS_I]);
+        e2i = full(sparse(range(0,max(i2e))', zeros(max(i2e), 1 ), zeros(max(i2e), 1 ), max(i2e), 1));
+        index_e2i = i2e - 1;
+        e2i = set(e2i, index_e2i, range(1, nb + 1 )' );
+
+        index_line = get_multi(e2i, slice(gen, [0], [GEN_BUS-1, GEN_BUS]) - 1) - 1;
+        Cld = full(sparse(index_line', range(0, ng)', is_ld, nb, ng));
+        if opt_nominal == 1{        // use nominal power
+            Pdd = -Cld * slice(gen, [0], [PMIN-1, PMIN]);      // real power
+
+            Q = zeros(ng, 1);
+            if ~~is_empty(ld){
+                value_Qld = (select(gen, [ld], [QMIN - 1]) == 0) .* select(gen, [ld], [QMAX - 1]) + (select(gen, [ld], [QMAX - 1]) == 0) .* select(gen, [ld], [QMIN - 1]);
+                Q = set(Q, ld, value_Qld );
+            }
+            Qdd = -Cld * Q;             // reactive power
+
+        }else{                // use realized actual power dispatch
+            Pdd = -Cld * slice(gen, [0], [PG-1, PG]);        // real power
+            Qdd = -Cld * slice(gen, [0], [QG-1, QG]);        // reactive power
+        }
+    }else{
+        Pdd = zeros(nb, 1);
+        Qdd = zeros(nb, 1);
+    }
+
+    // compute load sums
+    if (nz == nb) && (sum(LZone == range(1,nb+1)')/nz == 1){   //nz == nb且LZone和range(1,nb+1)'元素全相等
+        BUS_TYPE_NONE = (slice(bus, [0], [BUS_TYPE-1, BUS_TYPE]) == ones(nb,1)* NONE);  //判断节点类型是否为NONE
+        Pd = (Pdf + Pdd) .* (~~BUS_TYPE_NONE);
+        Qd = (Qdf + Qdd) .* (~~BUS_TYPE_NONE);
+    }else{
+        BUS_TYPE_NONE = (slice(bus, [0], [BUS_TYPE-1, BUS_TYPE]) == ones(nb,1)* NONE);
+        Pd = zeros(nz, 1);
+        Qd = zeros(nz, 1);
+
+        for k in 1..(nz+1){
+            idx = find( LZone == k && BUS_TYPE_NONE == 0);
+            Pd = set(Pd, [k-1], [sum(get_multi(Pdf, [idx])) + sum(get_multi(Pdd, [idx]))] );
+            Qd = set(Qd, [k-1], [sum(get_multi(Qdf, [idx])) + sum(get_multi(Qdd, [idx]))] );
+        }
+    }
 
     //return Pd, Qd;
     return Pd;