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