0001 function [f, df, d2f] = opf_costfcn(x, om, varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0039 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0040 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0041 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0042 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0043 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0044 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0045
0046
0047 mpc = get_mpc(om);
0048 [baseMVA, gen, gencost] = deal(mpc.baseMVA, mpc.gen, mpc.gencost);
0049 cp = get_cost_params(om);
0050 [N, Cw, H, dd, rh, kk, mm] = deal(cp.N, cp.Cw, cp.H, cp.dd, ...
0051 cp.rh, cp.kk, cp.mm);
0052 vv = get_idx(om);
0053
0054
0055 ng = size(gen, 1);
0056 ny = getN(om, 'var', 'y');
0057 nxyz = length(x);
0058
0059
0060 Pg = x(vv.i1.Pg:vv.iN.Pg);
0061 Qg = x(vv.i1.Qg:vv.iN.Qg);
0062
0063
0064
0065
0066
0067 ipol = find(gencost(:, MODEL) == POLYNOMIAL);
0068 xx = [ Pg; Qg ] * baseMVA;
0069 if ~isempty(ipol)
0070 f = sum( totcost(gencost(ipol, :), xx(ipol)) );
0071 else
0072 f = 0;
0073 end
0074
0075
0076 if ny > 0
0077 ccost = full(sparse(ones(1,ny), vv.i1.y:vv.iN.y, ones(1,ny), 1, nxyz));
0078 f = f + ccost * x;
0079 else
0080 ccost = zeros(1, nxyz);
0081 end
0082
0083
0084 if ~isempty(N)
0085 nw = size(N, 1);
0086 r = N * x - rh;
0087 iLT = find(r < -kk);
0088 iEQ = find(r == 0 & kk == 0);
0089 iGT = find(r > kk);
0090 iND = [iLT; iEQ; iGT];
0091 iL = find(dd == 1);
0092 iQ = find(dd == 2);
0093 LL = sparse(iL, iL, 1, nw, nw);
0094 QQ = sparse(iQ, iQ, 1, nw, nw);
0095 kbar = sparse(iND, iND, [ ones(length(iLT), 1);
0096 zeros(length(iEQ), 1);
0097 -ones(length(iGT), 1)], nw, nw) * kk;
0098 rr = r + kbar;
0099 M = sparse(iND, iND, mm(iND), nw, nw);
0100 diagrr = sparse(1:nw, 1:nw, rr, nw, nw);
0101
0102
0103 w = M * (LL + QQ * diagrr) * rr;
0104
0105 f = f + (w' * H * w) / 2 + Cw' * w;
0106 end
0107
0108
0109 if nargout > 1
0110
0111 iPg = vv.i1.Pg:vv.iN.Pg;
0112 iQg = vv.i1.Qg:vv.iN.Qg;
0113
0114
0115 df_dPgQg = zeros(2*ng, 1);
0116 df_dPgQg(ipol) = baseMVA * polycost(gencost(ipol, :), xx(ipol), 1);
0117 df = zeros(nxyz, 1);
0118 df(iPg) = df_dPgQg(1:ng);
0119 df(iQg) = df_dPgQg((1:ng) + ng);
0120
0121
0122 df = df + ccost';
0123
0124
0125
0126 if ~isempty(N)
0127 HwC = H * w + Cw;
0128 AA = N' * M * (LL + 2 * QQ * diagrr);
0129 df = df + AA * HwC;
0130
0131
0132 if 0
0133 ddff = zeros(size(df));
0134 step = 1e-7;
0135 tol = 1e-3;
0136 for k = 1:length(x)
0137 xx = x;
0138 xx(k) = xx(k) + step;
0139 ddff(k) = (opf_costfcn(xx, om) - f) / step;
0140 end
0141 if max(abs(ddff - df)) > tol
0142 idx = find(abs(ddff - df) == max(abs(ddff - df)));
0143 fprintf('\nMismatch in gradient\n');
0144 fprintf('idx df(num) df diff\n');
0145 fprintf('%4d%16g%16g%16g\n', [ 1:length(df); ddff'; df'; abs(ddff - df)' ]);
0146 fprintf('MAX\n');
0147 fprintf('%4d%16g%16g%16g\n', [ idx'; ddff(idx)'; df(idx)'; abs(ddff(idx) - df(idx))' ]);
0148 fprintf('\n');
0149 end
0150 end
0151 end
0152
0153
0154 if nargout > 2
0155 pcost = gencost(1:ng, :);
0156 if size(gencost, 1) > ng
0157 qcost = gencost(ng+1:2*ng, :);
0158 else
0159 qcost = [];
0160 end
0161
0162
0163 d2f_dPg2 = sparse(ng, 1);
0164 d2f_dQg2 = sparse(ng, 1);
0165 ipolp = find(pcost(:, MODEL) == POLYNOMIAL);
0166 d2f_dPg2(ipolp) = baseMVA^2 * polycost(pcost(ipolp, :), Pg(ipolp)*baseMVA, 2);
0167 if ~isempty(qcost)
0168 ipolq = find(qcost(:, MODEL) == POLYNOMIAL);
0169 d2f_dQg2(ipolq) = baseMVA^2 * polycost(qcost(ipolq, :), Qg(ipolq)*baseMVA, 2);
0170 end
0171 i = [iPg iQg]';
0172 d2f = sparse(i, i, [d2f_dPg2; d2f_dQg2], nxyz, nxyz);
0173
0174
0175 if ~isempty(N)
0176 d2f = d2f + AA * H * AA' + 2 * N' * M * QQ * sparse(1:nw, 1:nw, HwC, nw, nw) * N;
0177 end
0178 end
0179 end