Home > matpower5.1 > t > t_hessian.m

t_hessian

PURPOSE ^

T_HESSIAN Numerical tests of 2nd derivative code.

SYNOPSIS ^

function t_hessian(quiet)

DESCRIPTION ^

T_HESSIAN  Numerical tests of 2nd derivative code.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function t_hessian(quiet)
0002 %T_HESSIAN  Numerical tests of 2nd derivative code.
0003 
0004 %   MATPOWER
0005 %   Copyright (c) 2004-2015 by Power System Engineering Research Center (PSERC)
0006 %   by Ray Zimmerman, PSERC Cornell
0007 %
0008 %   $Id: t_hessian.m 2644 2015-03-11 19:34:22Z ray $
0009 %
0010 %   This file is part of MATPOWER.
0011 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0012 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0013 
0014 if nargin < 1
0015     quiet = 0;
0016 end
0017 
0018 t_begin(44, quiet);
0019 
0020 casefile = 'case30';
0021 
0022 %% define named indices into bus, gen, branch matrices
0023 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0024     VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0025 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0026     TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0027     ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0028 
0029 %% run powerflow to get solved case
0030 mpopt = mpoption('verbose', 0, 'out.all', 0);
0031 [baseMVA, bus, gen, branch, success, et] = runpf(casefile, mpopt);
0032 
0033 %% switch to internal bus numbering and build admittance matrices
0034 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0035 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0036 Vm = bus(:, VM);
0037 Va = bus(:, VA) * pi/180;
0038 V = Vm .* exp(1j * Va);
0039 f = branch(:, F_BUS);       %% list of "from" buses
0040 t = branch(:, T_BUS);       %% list of "to" buses
0041 nl = length(f);
0042 nb = length(V);
0043 Cf = sparse(1:nl, f, ones(nl, 1), nl, nb);      %% connection matrix for line & from buses
0044 Ct = sparse(1:nl, t, ones(nl, 1), nl, nb);      %% connection matrix for line & to buses
0045 pert = 1e-8;
0046 
0047 %%-----  check d2Sbus_dV2 code  -----
0048 t = ' - d2Sbus_dV2 (complex power injections)';
0049 lam = 10 * rand(nb, 1);
0050 num_Haa = zeros(nb, nb);
0051 num_Hav = zeros(nb, nb);
0052 num_Hva = zeros(nb, nb);
0053 num_Hvv = zeros(nb, nb);
0054 [dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
0055 [Haa, Hav, Hva, Hvv] = d2Sbus_dV2(Ybus, V, lam);
0056 for i = 1:nb
0057     Vap = V;
0058     Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0059     [dSbus_dVm_ap, dSbus_dVa_ap] = dSbus_dV(Ybus, Vap);
0060     num_Haa(:, i) = (dSbus_dVa_ap - dSbus_dVa).' * lam / pert;
0061     num_Hva(:, i) = (dSbus_dVm_ap - dSbus_dVm).' * lam / pert;
0062 
0063     Vmp = V;
0064     Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0065     [dSbus_dVm_mp, dSbus_dVa_mp] = dSbus_dV(Ybus, Vmp);
0066     num_Hav(:, i) = (dSbus_dVa_mp - dSbus_dVa).' * lam / pert;
0067     num_Hvv(:, i) = (dSbus_dVm_mp - dSbus_dVm).' * lam / pert;
0068 end
0069 
0070 t_is(full(Haa), num_Haa, 4, ['Haa' t]);
0071 t_is(full(Hav), num_Hav, 4, ['Hav' t]);
0072 t_is(full(Hva), num_Hva, 4, ['Hva' t]);
0073 t_is(full(Hvv), num_Hvv, 4, ['Hvv' t]);
0074 
0075 %%-----  check d2Sbr_dV2 code  -----
0076 t = ' - d2Sbr_dV2 (complex power flows)';
0077 lam = 10 * rand(nl, 1);
0078 % lam = [1; zeros(nl-1, 1)];
0079 num_Gfaa = zeros(nb, nb);
0080 num_Gfav = zeros(nb, nb);
0081 num_Gfva = zeros(nb, nb);
0082 num_Gfvv = zeros(nb, nb);
0083 num_Gtaa = zeros(nb, nb);
0084 num_Gtav = zeros(nb, nb);
0085 num_Gtva = zeros(nb, nb);
0086 num_Gtvv = zeros(nb, nb);
0087 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0088 [Gfaa, Gfav, Gfva, Gfvv] = d2Sbr_dV2(Cf, Yf, V, lam);
0089 [Gtaa, Gtav, Gtva, Gtvv] = d2Sbr_dV2(Ct, Yt, V, lam);
0090 for i = 1:nb
0091     Vap = V;
0092     Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0093     [dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap] = ...
0094         dSbr_dV(branch, Yf, Yt, Vap);
0095     num_Gfaa(:, i) = (dSf_dVa_ap - dSf_dVa).' * lam / pert;
0096     num_Gfva(:, i) = (dSf_dVm_ap - dSf_dVm).' * lam / pert;
0097     num_Gtaa(:, i) = (dSt_dVa_ap - dSt_dVa).' * lam / pert;
0098     num_Gtva(:, i) = (dSt_dVm_ap - dSt_dVm).' * lam / pert;
0099 
0100     Vmp = V;
0101     Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0102     [dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp] = ...
0103         dSbr_dV(branch, Yf, Yt, Vmp);
0104     num_Gfav(:, i) = (dSf_dVa_mp - dSf_dVa).' * lam / pert;
0105     num_Gfvv(:, i) = (dSf_dVm_mp - dSf_dVm).' * lam / pert;
0106     num_Gtav(:, i) = (dSt_dVa_mp - dSt_dVa).' * lam / pert;
0107     num_Gtvv(:, i) = (dSt_dVm_mp - dSt_dVm).' * lam / pert;
0108 end
0109 
0110 t_is(full(Gfaa), num_Gfaa, 4, ['Gfaa' t]);
0111 t_is(full(Gfav), num_Gfav, 4, ['Gfav' t]);
0112 t_is(full(Gfva), num_Gfva, 4, ['Gfva' t]);
0113 t_is(full(Gfvv), num_Gfvv, 4, ['Gfvv' t]);
0114 
0115 t_is(full(Gtaa), num_Gtaa, 4, ['Gtaa' t]);
0116 t_is(full(Gtav), num_Gtav, 4, ['Gtav' t]);
0117 t_is(full(Gtva), num_Gtva, 4, ['Gtva' t]);
0118 t_is(full(Gtvv), num_Gtvv, 4, ['Gtvv' t]);
0119 
0120 %%-----  check d2Ibr_dV2 code  -----
0121 t = ' - d2Ibr_dV2 (complex currents)';
0122 lam = 10 * rand(nl, 1);
0123 % lam = [1; zeros(nl-1, 1)];
0124 num_Gfaa = zeros(nb, nb);
0125 num_Gfav = zeros(nb, nb);
0126 num_Gfva = zeros(nb, nb);
0127 num_Gfvv = zeros(nb, nb);
0128 num_Gtaa = zeros(nb, nb);
0129 num_Gtav = zeros(nb, nb);
0130 num_Gtva = zeros(nb, nb);
0131 num_Gtvv = zeros(nb, nb);
0132 [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch, Yf, Yt, V);
0133 [Gfaa, Gfav, Gfva, Gfvv] = d2Ibr_dV2(Yf, V, lam);
0134 [Gtaa, Gtav, Gtva, Gtvv] = d2Ibr_dV2(Yt, V, lam);
0135 for i = 1:nb
0136     Vap = V;
0137     Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0138     [dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap] = ...
0139         dIbr_dV(branch, Yf, Yt, Vap);
0140     num_Gfaa(:, i) = (dIf_dVa_ap - dIf_dVa).' * lam / pert;
0141     num_Gfva(:, i) = (dIf_dVm_ap - dIf_dVm).' * lam / pert;
0142     num_Gtaa(:, i) = (dIt_dVa_ap - dIt_dVa).' * lam / pert;
0143     num_Gtva(:, i) = (dIt_dVm_ap - dIt_dVm).' * lam / pert;
0144 
0145     Vmp = V;
0146     Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0147     [dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp] = ...
0148         dIbr_dV(branch, Yf, Yt, Vmp);
0149     num_Gfav(:, i) = (dIf_dVa_mp - dIf_dVa).' * lam / pert;
0150     num_Gfvv(:, i) = (dIf_dVm_mp - dIf_dVm).' * lam / pert;
0151     num_Gtav(:, i) = (dIt_dVa_mp - dIt_dVa).' * lam / pert;
0152     num_Gtvv(:, i) = (dIt_dVm_mp - dIt_dVm).' * lam / pert;
0153 end
0154 
0155 t_is(full(Gfaa), num_Gfaa, 4, ['Gfaa' t]);
0156 t_is(full(Gfav), num_Gfav, 4, ['Gfav' t]);
0157 t_is(full(Gfva), num_Gfva, 4, ['Gfva' t]);
0158 t_is(full(Gfvv), num_Gfvv, 4, ['Gfvv' t]);
0159 
0160 t_is(full(Gtaa), num_Gtaa, 4, ['Gtaa' t]);
0161 t_is(full(Gtav), num_Gtav, 4, ['Gtav' t]);
0162 t_is(full(Gtva), num_Gtva, 4, ['Gtva' t]);
0163 t_is(full(Gtvv), num_Gtvv, 4, ['Gtvv' t]);
0164 
0165 %%-----  check d2ASbr_dV2 code  -----
0166 t = ' - d2ASbr_dV2 (squared apparent power flows)';
0167 lam = 10 * rand(nl, 1);
0168 % lam = [1; zeros(nl-1, 1)];
0169 num_Gfaa = zeros(nb, nb);
0170 num_Gfav = zeros(nb, nb);
0171 num_Gfva = zeros(nb, nb);
0172 num_Gfvv = zeros(nb, nb);
0173 num_Gtaa = zeros(nb, nb);
0174 num_Gtav = zeros(nb, nb);
0175 num_Gtva = zeros(nb, nb);
0176 num_Gtvv = zeros(nb, nb);
0177 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0178 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0179                         dAbr_dV(dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St);
0180 [Gfaa, Gfav, Gfva, Gfvv] = d2ASbr_dV2(dSf_dVa, dSf_dVm, Sf, Cf, Yf, V, lam);
0181 [Gtaa, Gtav, Gtva, Gtvv] = d2ASbr_dV2(dSt_dVa, dSt_dVm, St, Ct, Yt, V, lam);
0182 for i = 1:nb
0183     Vap = V;
0184     Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0185     [dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap] = ...
0186         dSbr_dV(branch, Yf, Yt, Vap);
0187     [dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap] = ...
0188         dAbr_dV(dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap);
0189     num_Gfaa(:, i) = (dAf_dVa_ap - dAf_dVa).' * lam / pert;
0190     num_Gfva(:, i) = (dAf_dVm_ap - dAf_dVm).' * lam / pert;
0191     num_Gtaa(:, i) = (dAt_dVa_ap - dAt_dVa).' * lam / pert;
0192     num_Gtva(:, i) = (dAt_dVm_ap - dAt_dVm).' * lam / pert;
0193 
0194     Vmp = V;
0195     Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0196     [dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp] = ...
0197         dSbr_dV(branch, Yf, Yt, Vmp);
0198     [dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp] = ...
0199         dAbr_dV(dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp);
0200     num_Gfav(:, i) = (dAf_dVa_mp - dAf_dVa).' * lam / pert;
0201     num_Gfvv(:, i) = (dAf_dVm_mp - dAf_dVm).' * lam / pert;
0202     num_Gtav(:, i) = (dAt_dVa_mp - dAt_dVa).' * lam / pert;
0203     num_Gtvv(:, i) = (dAt_dVm_mp - dAt_dVm).' * lam / pert;
0204 end
0205 
0206 t_is(full(Gfaa), num_Gfaa, 2, ['Gfaa' t]);
0207 t_is(full(Gfav), num_Gfav, 2, ['Gfav' t]);
0208 t_is(full(Gfva), num_Gfva, 2, ['Gfva' t]);
0209 t_is(full(Gfvv), num_Gfvv, 2, ['Gfvv' t]);
0210 
0211 t_is(full(Gtaa), num_Gtaa, 2, ['Gtaa' t]);
0212 t_is(full(Gtav), num_Gtav, 2, ['Gtav' t]);
0213 t_is(full(Gtva), num_Gtva, 2, ['Gtva' t]);
0214 t_is(full(Gtvv), num_Gtvv, 2, ['Gtvv' t]);
0215 
0216 %%-----  check d2ASbr_dV2 code  -----
0217 t = ' - d2ASbr_dV2 (squared real power flows)';
0218 lam = 10 * rand(nl, 1);
0219 % lam = [1; zeros(nl-1, 1)];
0220 num_Gfaa = zeros(nb, nb);
0221 num_Gfav = zeros(nb, nb);
0222 num_Gfva = zeros(nb, nb);
0223 num_Gfvv = zeros(nb, nb);
0224 num_Gtaa = zeros(nb, nb);
0225 num_Gtav = zeros(nb, nb);
0226 num_Gtva = zeros(nb, nb);
0227 num_Gtvv = zeros(nb, nb);
0228 [dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
0229 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0230                         dAbr_dV(real(dSf_dVa), real(dSf_dVm), real(dSt_dVa), real(dSt_dVm), real(Sf), real(St));
0231 [Gfaa, Gfav, Gfva, Gfvv] = d2ASbr_dV2(real(dSf_dVa), real(dSf_dVm), real(Sf), Cf, Yf, V, lam);
0232 [Gtaa, Gtav, Gtva, Gtvv] = d2ASbr_dV2(real(dSt_dVa), real(dSt_dVm), real(St), Ct, Yt, V, lam);
0233 for i = 1:nb
0234     Vap = V;
0235     Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0236     [dSf_dVa_ap, dSf_dVm_ap, dSt_dVa_ap, dSt_dVm_ap, Sf_ap, St_ap] = ...
0237         dSbr_dV(branch, Yf, Yt, Vap);
0238     [dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap] = ...
0239         dAbr_dV(real(dSf_dVa_ap), real(dSf_dVm_ap), real(dSt_dVa_ap), real(dSt_dVm_ap), real(Sf_ap), real(St_ap));
0240     num_Gfaa(:, i) = (dAf_dVa_ap - dAf_dVa).' * lam / pert;
0241     num_Gfva(:, i) = (dAf_dVm_ap - dAf_dVm).' * lam / pert;
0242     num_Gtaa(:, i) = (dAt_dVa_ap - dAt_dVa).' * lam / pert;
0243     num_Gtva(:, i) = (dAt_dVm_ap - dAt_dVm).' * lam / pert;
0244 
0245     Vmp = V;
0246     Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0247     [dSf_dVa_mp, dSf_dVm_mp, dSt_dVa_mp, dSt_dVm_mp, Sf_mp, St_mp] = ...
0248         dSbr_dV(branch, Yf, Yt, Vmp);
0249     [dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp] = ...
0250         dAbr_dV(real(dSf_dVa_mp), real(dSf_dVm_mp), real(dSt_dVa_mp), real(dSt_dVm_mp), real(Sf_mp), real(St_mp));
0251     num_Gfav(:, i) = (dAf_dVa_mp - dAf_dVa).' * lam / pert;
0252     num_Gfvv(:, i) = (dAf_dVm_mp - dAf_dVm).' * lam / pert;
0253     num_Gtav(:, i) = (dAt_dVa_mp - dAt_dVa).' * lam / pert;
0254     num_Gtvv(:, i) = (dAt_dVm_mp - dAt_dVm).' * lam / pert;
0255 end
0256 
0257 t_is(full(Gfaa), num_Gfaa, 2, ['Gfaa' t]);
0258 t_is(full(Gfav), num_Gfav, 2, ['Gfav' t]);
0259 t_is(full(Gfva), num_Gfva, 2, ['Gfva' t]);
0260 t_is(full(Gfvv), num_Gfvv, 2, ['Gfvv' t]);
0261 
0262 t_is(full(Gtaa), num_Gtaa, 2, ['Gtaa' t]);
0263 t_is(full(Gtav), num_Gtav, 2, ['Gtav' t]);
0264 t_is(full(Gtva), num_Gtva, 2, ['Gtva' t]);
0265 t_is(full(Gtvv), num_Gtvv, 2, ['Gtvv' t]);
0266 
0267 %%-----  check d2AIbr_dV2 code  -----
0268 t = ' - d2AIbr_dV2 (squared current magnitudes)';
0269 lam = 10 * rand(nl, 1);
0270 % lam = [1; zeros(nl-1, 1)];
0271 num_Gfaa = zeros(nb, nb);
0272 num_Gfav = zeros(nb, nb);
0273 num_Gfva = zeros(nb, nb);
0274 num_Gfvv = zeros(nb, nb);
0275 num_Gtaa = zeros(nb, nb);
0276 num_Gtav = zeros(nb, nb);
0277 num_Gtva = zeros(nb, nb);
0278 num_Gtvv = zeros(nb, nb);
0279 [dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It] = dIbr_dV(branch, Yf, Yt, V);
0280 [dAf_dVa, dAf_dVm, dAt_dVa, dAt_dVm] = ...
0281                         dAbr_dV(dIf_dVa, dIf_dVm, dIt_dVa, dIt_dVm, If, It);
0282 [Gfaa, Gfav, Gfva, Gfvv] = d2AIbr_dV2(dIf_dVa, dIf_dVm, If, Yf, V, lam);
0283 [Gtaa, Gtav, Gtva, Gtvv] = d2AIbr_dV2(dIt_dVa, dIt_dVm, It, Yt, V, lam);
0284 for i = 1:nb
0285     Vap = V;
0286     Vap(i) = Vm(i) * exp(1j * (Va(i) + pert));
0287     [dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap] = ...
0288         dIbr_dV(branch, Yf, Yt, Vap);
0289     [dAf_dVa_ap, dAf_dVm_ap, dAt_dVa_ap, dAt_dVm_ap] = ...
0290         dAbr_dV(dIf_dVa_ap, dIf_dVm_ap, dIt_dVa_ap, dIt_dVm_ap, If_ap, It_ap);
0291     num_Gfaa(:, i) = (dAf_dVa_ap - dAf_dVa).' * lam / pert;
0292     num_Gfva(:, i) = (dAf_dVm_ap - dAf_dVm).' * lam / pert;
0293     num_Gtaa(:, i) = (dAt_dVa_ap - dAt_dVa).' * lam / pert;
0294     num_Gtva(:, i) = (dAt_dVm_ap - dAt_dVm).' * lam / pert;
0295 
0296     Vmp = V;
0297     Vmp(i) = (Vm(i) + pert) * exp(1j * Va(i));
0298     [dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp] = ...
0299         dIbr_dV(branch, Yf, Yt, Vmp);
0300     [dAf_dVa_mp, dAf_dVm_mp, dAt_dVa_mp, dAt_dVm_mp] = ...
0301         dAbr_dV(dIf_dVa_mp, dIf_dVm_mp, dIt_dVa_mp, dIt_dVm_mp, If_mp, It_mp);
0302     num_Gfav(:, i) = (dAf_dVa_mp - dAf_dVa).' * lam / pert;
0303     num_Gfvv(:, i) = (dAf_dVm_mp - dAf_dVm).' * lam / pert;
0304     num_Gtav(:, i) = (dAt_dVa_mp - dAt_dVa).' * lam / pert;
0305     num_Gtvv(:, i) = (dAt_dVm_mp - dAt_dVm).' * lam / pert;
0306 end
0307 
0308 t_is(full(Gfaa), num_Gfaa, 3, ['Gfaa' t]);
0309 t_is(full(Gfav), num_Gfav, 3, ['Gfav' t]);
0310 t_is(full(Gfva), num_Gfva, 3, ['Gfva' t]);
0311 t_is(full(Gfvv), num_Gfvv, 2, ['Gfvv' t]);
0312 
0313 t_is(full(Gtaa), num_Gtaa, 3, ['Gtaa' t]);
0314 t_is(full(Gtav), num_Gtav, 3, ['Gtav' t]);
0315 t_is(full(Gtva), num_Gtva, 3, ['Gtva' t]);
0316 t_is(full(Gtvv), num_Gtvv, 2, ['Gtvv' t]);
0317 
0318 t_end;
Generated on Fri 20-Mar-2015 18:23:34 by m2html © 2005