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