0001 function om = opf_setup(mpc, mpopt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 dc = strcmp(upper(mpopt.model), 'DC');
0021 alg = upper(mpopt.opf.ac.solver);
0022 use_vg = mpopt.opf.use_vg;
0023 vcart = ~dc && mpopt.opf.v_cartesian;
0024
0025
0026 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0027 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0028 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0029 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0030 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0031 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0032 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0033 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0034 [PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
0035
0036
0037
0038
0039 if strcmp(alg, 'MINOPF') || strcmp(alg, 'PDIPM') || ...
0040 strcmp(alg, 'TRALM') || strcmp(alg, 'SDPOPF')
0041 legacy_formulation = 1;
0042 if vcart
0043 error('Option ''opf.v_cartesian'' = 1 is not compatible with ''opf.solver.ac''=''%s''.', alg);
0044 end
0045 if mpopt.opf.current_balance
0046 error('Option ''opf.current_balance'' = 1 is not compatible with ''opf.solver.ac''=''%s''.', alg);
0047 end
0048 else
0049 legacy_formulation = 0;
0050 end
0051 if ~dc && ( ~isempty(mpopt.exp.sys_wide_zip_loads.pw) && ...
0052 ~isequal(mpopt.exp.sys_wide_zip_loads.pw, [1 0 0]) || ...
0053 ~isempty(mpopt.exp.sys_wide_zip_loads.qw) && ...
0054 ~isequal(mpopt.exp.sys_wide_zip_loads.qw, [1 0 0]) )
0055 if vcart
0056 warning('Voltage dependent loads are not supported with option ''opf.v_cartesian'' = 1. Reverting to constant power load model.');
0057 end
0058 if mpopt.opf.current_balance
0059 warning('Voltage dependent loads are not supported with option ''opf.current_balance'' = 1. Reverting to constant power load model.');
0060 end
0061 end
0062
0063
0064 nb = size(mpc.bus, 1);
0065 nl = size(mpc.branch, 1);
0066 ng = size(mpc.gen, 1);
0067 nnle = 0;
0068 nnli = 0;
0069 if isfield(mpc, 'A')
0070 nlin = size(mpc.A, 1);
0071 else
0072 nlin = 0;
0073 end
0074 if isfield(mpc, 'N')
0075 nw = size(mpc.N, 1);
0076 else
0077 nw = 0;
0078 end
0079
0080 if dc
0081
0082 mpc.gencost = pqcost(mpc.gencost, ng);
0083
0084
0085 if nlin || nw
0086 acc = [nb+(1:nb) 2*nb+ng+(1:ng)];
0087 if nlin && size(mpc.A, 2) >= 2*nb + 2*ng
0088
0089 if any(any(mpc.A(:, acc)))
0090 error('opf_setup: attempting to solve DC OPF with user constraints on Vm or Qg');
0091 end
0092 mpc.A(:, acc) = [];
0093 end
0094 if nw && size(mpc.N, 2) >= 2*nb + 2*ng
0095
0096 if any(any(mpc.N(:, acc)))
0097 [ii, jj] = find(mpc.N(:, acc));
0098 ii = unique(ii);
0099 if any(mpc.Cw(ii)) || (isfield(mpc, 'H') && ~isempty(mpc.H) && ...
0100 any(any(mpc.H(:, ii))))
0101 error('opf_setup: attempting to solve DC OPF with user costs on Vm or Qg');
0102 end
0103 end
0104 mpc.N(:, acc) = [];
0105 end
0106 end
0107 else
0108 if use_vg
0109
0110 Cg = sparse(mpc.gen(:, GEN_BUS), (1:ng)', mpc.gen(:, GEN_STATUS) > 0, nb, ng);
0111 Vbg = Cg * sparse(1:ng, 1:ng, mpc.gen(:, VG), ng, ng);
0112 Vmax = max(Vbg, [], 2);
0113 ib = find(Vmax);
0114 Vmin = max(2*Cg - Vbg, [], 2);
0115 Vmin(ib) = 2 - Vmin(ib);
0116
0117 if use_vg == 1
0118 mpc.bus(ib, VMAX) = Vmax(ib);
0119 mpc.bus(ib, VMIN) = Vmin(ib);
0120 mpc.bus(ib, VM) = mpc.bus(ib, VMAX);
0121 elseif use_vg > 0 && use_vg < 1
0122
0123 mpc.bus(ib, VMAX) = (1-use_vg) * mpc.bus(ib, VMAX) + use_vg * Vmax(ib);
0124 mpc.bus(ib, VMIN) = (1-use_vg) * mpc.bus(ib, VMIN) + use_vg * Vmin(ib);
0125 else
0126 error('opf_setup: option ''opf.use_vg'' (= %g) cannot be negative or greater than 1', use_vg);
0127 end
0128 end
0129 if isfield(mpc, 'user_constraints')
0130 if isfield(mpc.user_constraints, 'nle')
0131 for k = 1:length(mpc.user_constraints.nle)
0132 nnle = nnle + mpc.user_constraints.nle{k}{2};
0133 end
0134 end
0135 if isfield(mpc.user_constraints, 'nli')
0136 for k = 1:length(mpc.user_constraints.nli)
0137 nnli = nnli + mpc.user_constraints.nli{k}{2};
0138 end
0139 end
0140 end
0141 end
0142
0143
0144 pwl1 = find(mpc.gencost(:, MODEL) == PW_LINEAR & mpc.gencost(:, NCOST) == 2);
0145
0146 if ~isempty(pwl1)
0147 x0 = mpc.gencost(pwl1, COST);
0148 y0 = mpc.gencost(pwl1, COST+1);
0149 x1 = mpc.gencost(pwl1, COST+2);
0150 y1 = mpc.gencost(pwl1, COST+3);
0151 m = (y1 - y0) ./ (x1 - x0);
0152 b = y0 - m .* x0;
0153 mpc.gencost(pwl1, MODEL) = POLYNOMIAL;
0154 mpc.gencost(pwl1, NCOST) = 2;
0155 mpc.gencost(pwl1, COST:COST+1) = [m b];
0156 end
0157
0158
0159 [baseMVA, bus, gen, branch, gencost, Au, lbu, ubu, mpopt, ...
0160 N, fparm, H, Cw, z0, zl, zu, userfcn] = opf_args(mpc, mpopt);
0161
0162
0163 refs = find(bus(:, BUS_TYPE) == REF);
0164 if length(refs) > 1 && mpopt.verbose > 0
0165 errstr = ['\nopf_setup: Warning: Multiple reference buses.\n', ...
0166 ' For a system with islands, a reference bus in each island\n', ...
0167 ' may help convergence, but in a fully connected system such\n', ...
0168 ' a situation is probably not reasonable.\n\n' ];
0169 fprintf(errstr);
0170 end
0171
0172
0173 Va = bus(:, VA) * (pi/180);
0174 Vau = Inf(nb, 1);
0175 Val = -Vau;
0176 Vau(refs) = Va(refs);
0177 Val(refs) = Va(refs);
0178 Pg = gen(:, PG) / baseMVA;
0179 Pmin = gen(:, PMIN) / baseMVA;
0180 Pmax = gen(:, PMAX) / baseMVA;
0181 if ~dc
0182 Vm = bus(:, VM);
0183 Qg = gen(:, QG) / baseMVA;
0184 Qmin = gen(:, QMIN) / baseMVA;
0185 Qmax = gen(:, QMAX) / baseMVA;
0186 if vcart
0187 V = Vm .* exp(1j*Va);
0188 Vr = real(V);
0189 Vi = imag(V);
0190 end
0191 end
0192
0193
0194 cpg = [];
0195 cqg = [];
0196 [pcost qcost] = pqcost(mpc.gencost, ng);
0197 ip0 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 1);
0198 ip1 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 2);
0199 ip2 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) == 3);
0200 ip3 = find(pcost(:, MODEL) == POLYNOMIAL & pcost(:, NCOST) > 3);
0201 if ~isempty(ip2) || ~isempty(ip1) || ~isempty(ip0)
0202 kpg = zeros(ng, 1);
0203 cpg = zeros(ng, 1);
0204 if ~isempty(ip2)
0205 Qpg = zeros(ng, 1);
0206 Qpg(ip2) = 2 * pcost(ip2, COST) * baseMVA^2;
0207 cpg(ip2) = cpg(ip2) + pcost(ip2, COST+1) * baseMVA;
0208 kpg(ip2) = kpg(ip2) + pcost(ip2, COST+2);
0209 else
0210 Qpg = [];
0211 end
0212 if ~isempty(ip1)
0213 cpg(ip1) = cpg(ip1) + pcost(ip1, COST) * baseMVA;
0214 kpg(ip1) = kpg(ip1) + pcost(ip1, COST+1);
0215 end
0216 if ~isempty(ip0)
0217 kpg(ip0) = kpg(ip0) + pcost(ip0, COST);
0218 end
0219 end
0220 if ~isempty(qcost)
0221 iq0 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 1);
0222 iq1 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 2);
0223 iq2 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) == 3);
0224 iq3 = find(qcost(:, MODEL) == POLYNOMIAL & qcost(:, NCOST) > 3);
0225 if ~isempty(iq2) || ~isempty(iq1) || ~isempty(iq0)
0226 kqg = zeros(ng, 1);
0227 cqg = zeros(ng, 1);
0228 if ~isempty(iq2)
0229 Qqg = zeros(ng, 1);
0230 Qqg(iq2) = 2 * qcost(iq2, COST) * baseMVA^2;
0231 cqg(iq2) = cqg(iq2) + qcost(iq2, COST+1) * baseMVA;
0232 kqg(iq2) = kqg(iq2) + qcost(iq2, COST+2);
0233 else
0234 Qqg = [];
0235 end
0236 if ~isempty(iq1)
0237 cqg(iq1) = cqg(iq1) + qcost(iq1, COST) * baseMVA;
0238 kqg(iq1) = kqg(iq1) + qcost(iq1, COST+1);
0239 end
0240 if ~isempty(iq0)
0241 kqg(iq0) = kqg(iq0) + qcost(iq0, COST);
0242 end
0243 end
0244 end
0245
0246
0247 [Aang, lang, uang, iang] = makeAang(baseMVA, branch, nb, mpopt);
0248 nang = length(iang);
0249
0250 if dc
0251
0252 if ~isempty(ip3)
0253 error('opf_setup: DC OPF cannot handle polynomial costs with higher than quadratic order.');
0254 end
0255
0256
0257 nv = 0;
0258 nq = 0;
0259 q1 = [];
0260
0261
0262 [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0263 neg_Cg = sparse(gen(:, GEN_BUS), 1:ng, -1, nb, ng);
0264 Amis = [B neg_Cg];
0265 bmis = -(bus(:, PD) + bus(:, GS)) / baseMVA - Pbusinj;
0266
0267
0268 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0269 nl2 = length(il);
0270 if nl2
0271 upf = branch(il, RATE_A) / baseMVA - Pfinj(il);
0272 upt = branch(il, RATE_A) / baseMVA + Pfinj(il);
0273 else
0274 upf = [];
0275 upt = [];
0276 end
0277
0278 user_vars = {'Va', 'Pg'};
0279 ycon_vars = {'Pg', 'y'};
0280 else
0281
0282 nv = nb;
0283 nq = ng;
0284 q1 = 1+ng;
0285
0286
0287 il = find(branch(:, RATE_A) ~= 0 & branch(:, RATE_A) < 1e10);
0288 nl2 = length(il);
0289
0290
0291 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0292
0293
0294 [Avl, lvl, uvl] = makeAvl(mpc);
0295
0296
0297 [Apqh, ubpqh, Apql, ubpql, Apqdata] = makeApq(baseMVA, gen);
0298
0299 if vcart
0300 user_vars = {'Vr', 'Vi', 'Pg', 'Qg'};
0301 nodal_balance_vars = {'Vr', 'Vi', 'Pg', 'Qg'};
0302 flow_lim_vars = {'Vr', 'Vi'};
0303 else
0304 user_vars = {'Va', 'Vm', 'Pg', 'Qg'};
0305 nodal_balance_vars = {'Va', 'Vm', 'Pg', 'Qg'};
0306 flow_lim_vars = {'Va', 'Vm'};
0307 end
0308 ycon_vars = {'Pg', 'Qg', 'y'};
0309
0310
0311 if mpopt.opf.current_balance
0312 mis_cons = {'rImis', 'iImis'};
0313 fcn_mis = @(x)opf_current_balance_fcn(x, mpc, Ybus, mpopt);
0314 hess_mis = @(x, lam)opf_current_balance_hess(x, lam, mpc, Ybus, mpopt);
0315 else
0316 mis_cons = {'Pmis', 'Qmis'};
0317 fcn_mis = @(x)opf_power_balance_fcn(x, mpc, Ybus, mpopt);
0318 hess_mis = @(x, lam)opf_power_balance_hess(x, lam, mpc, Ybus, mpopt);
0319 end
0320 fcn_flow = @(x)opf_branch_flow_fcn(x, mpc, Yf(il, :), Yt(il, :), il, mpopt);
0321 hess_flow = @(x, lam)opf_branch_flow_hess(x, lam, mpc, Yf(il, :), Yt(il, :), il, mpopt);
0322 if vcart
0323 fcn_vref = @(x)opf_vref_fcn(x, mpc, refs, mpopt);
0324 hess_vref = @(x, lam)opf_vref_hess(x, lam, mpc, refs, mpopt);
0325 veq = find(mpc.bus(:, VMIN) == mpc.bus(:, VMAX));
0326 viq = find(mpc.bus(:, VMIN) ~= mpc.bus(:, VMAX));
0327 nveq = length(veq);
0328 nvlims = length(viq);
0329 if nveq
0330 fcn_veq = @(x)opf_veq_fcn(x, mpc, veq, mpopt);
0331 hess_veq = @(x, lam)opf_veq_hess(x, lam, mpc, veq, mpopt);
0332 end
0333 fcn_vlim = @(x)opf_vlim_fcn(x, mpc, viq, mpopt);
0334 hess_vlim = @(x, lam)opf_vlim_hess(x, lam, mpc, viq, mpopt);
0335 fcn_ang = @(x)opf_branch_ang_fcn(x, Aang, lang, uang);
0336 hess_ang = @(x, lam)opf_branch_ang_hess(x, lam, Aang, lang, uang);
0337 end
0338
0339
0340 if ~isempty(ip3)
0341 cost_Pg = @(x)opf_gen_cost_fcn(x, baseMVA, pcost, ip3, mpopt);
0342 end
0343 if ~isempty(qcost) && ~isempty(iq3)
0344 cost_Qg = @(x)opf_gen_cost_fcn(x, baseMVA, qcost, iq3, mpopt);
0345 end
0346 end
0347
0348
0349 if (strcmp(alg, 'PDIPM') && mpopt.pdipm.step_control) || strcmp(alg, 'TRALM')
0350
0351 ny = 0;
0352 Ay = sparse(0, ng+nq);
0353 by =[];
0354 else
0355 ipwl = find(gencost(:, MODEL) == PW_LINEAR);
0356 ny = size(ipwl, 1);
0357 [Ay, by] = makeAy(baseMVA, ng, gencost, 1, q1, 1+ng+nq);
0358 end
0359 if any(gencost(:, MODEL) ~= POLYNOMIAL & gencost(:, MODEL) ~= PW_LINEAR)
0360 error('opf_setup: some generator cost rows have invalid MODEL value');
0361 end
0362
0363
0364 nx = nb+nv + ng+nq;
0365 if nlin
0366 nz = size(mpc.A, 2) - nx;
0367 if nz < 0
0368 error('opf_setup: user supplied A matrix must have at least %d columns.', nx);
0369 end
0370 else
0371 nz = 0;
0372 if nw
0373 if size(mpc.N, 2) ~= nx;
0374 error('opf_setup: user supplied N matrix must have %d columns.', nx);
0375 end
0376 end
0377 end
0378
0379
0380 om = opf_model(mpc);
0381 if ~isempty(pwl1)
0382 om.userdata.pwl1 = pwl1;
0383 end
0384 om.userdata.iang = iang;
0385 if dc
0386
0387 om.userdata.Bf = Bf;
0388 om.userdata.Pfinj = Pfinj;
0389
0390
0391 om.add_var('Va', nb, Va, Val, Vau);
0392 om.add_var('Pg', ng, Pg, Pmin, Pmax);
0393
0394
0395 om.add_lin_constraint('Pmis', Amis, bmis, bmis, {'Va', 'Pg'});
0396 om.add_lin_constraint('Pf', Bf(il,:), -upt, upf, {'Va'});
0397 om.add_lin_constraint('ang', Aang, lang, uang, {'Va'});
0398
0399
0400 if ~isempty(cpg)
0401 om.add_quad_cost('polPg', Qpg, cpg, kpg, {'Pg'});
0402 end
0403 else
0404
0405 om.userdata.Apqdata = Apqdata;
0406
0407
0408 if vcart
0409 Vclim = 1.1 * bus(:, VMAX);
0410 om.add_var('Vr', nb, Vr, -Vclim, Vclim);
0411 om.add_var('Vi', nb, Vi, -Vclim, Vclim);
0412 else
0413 om.add_var('Va', nb, Va, Val, Vau);
0414 om.add_var('Vm', nb, Vm, bus(:, VMIN), bus(:, VMAX));
0415 end
0416 om.add_var('Pg', ng, Pg, Pmin, Pmax);
0417 om.add_var('Qg', ng, Qg, Qmin, Qmax);
0418
0419
0420 om.add_nln_constraint(mis_cons, [nb;nb], 1, fcn_mis, hess_mis, nodal_balance_vars);
0421 if legacy_formulation
0422 om.add_nln_constraint({'Sf', 'St'}, [nl;nl], 0, fcn_flow, hess_flow, flow_lim_vars);
0423 else
0424 om.add_nln_constraint({'Sf', 'St'}, [nl2;nl2], 0, fcn_flow, hess_flow, flow_lim_vars);
0425 end
0426 if vcart
0427 om.userdata.veq = veq;
0428 om.userdata.viq = viq;
0429 om.add_nln_constraint('Vref', length(refs), 1, fcn_vref, hess_vref, {'Vr', 'Vi'});
0430 if nveq
0431 om.add_nln_constraint('Veq', nveq, 1, fcn_veq, hess_veq, {'Vr', 'Vi'});
0432 end
0433 om.add_nln_constraint({'Vmin', 'Vmax'}, [nvlims;nvlims], 0, fcn_vlim, hess_vlim, {'Vr', 'Vi'});
0434 om.add_nln_constraint({'angL', 'angU'}, [nang;nang], 0, fcn_ang, hess_ang, {'Vr', 'Vi'});
0435 end
0436
0437
0438 om.add_lin_constraint('PQh', Apqh, [], ubpqh, {'Pg', 'Qg'});
0439 om.add_lin_constraint('PQl', Apql, [], ubpql, {'Pg', 'Qg'});
0440 om.add_lin_constraint('vl', Avl, lvl, uvl, {'Pg', 'Qg'});
0441 if ~vcart
0442 om.add_lin_constraint('ang', Aang, lang, uang, {'Va'});
0443 end
0444
0445
0446 if ~legacy_formulation
0447
0448 if ~isempty(cpg)
0449 om.add_quad_cost('polPg', Qpg, cpg, kpg, {'Pg'});
0450 end
0451 if ~isempty(cqg)
0452 om.add_quad_cost('polQg', Qqg, cqg, kqg, {'Qg'});
0453 end
0454
0455
0456 if ~isempty(ip3)
0457 om.add_nln_cost('polPg', 1, cost_Pg, {'Pg'});
0458 end
0459 if ~isempty(qcost) && ~isempty(iq3)
0460 om.add_nln_cost('polQg', 1, cost_Qg, {'Qg'});
0461 end
0462 end
0463 end
0464
0465
0466 if ny > 0
0467 om.add_var('y', ny);
0468 om.add_lin_constraint('ycon', Ay, [], by, ycon_vars);
0469 if dc || ~legacy_formulation
0470 om.add_quad_cost('pwl', [], ones(ny, 1), 0, {'y'});
0471 end
0472 end
0473
0474
0475 if nz > 0
0476 om.add_var('z', nz, z0, zl, zu);
0477 user_vars{end+1} = 'z';
0478 end
0479 if nlin
0480 om.add_lin_constraint('usr', mpc.A, lbu, ubu, user_vars);
0481 end
0482 if nnle
0483 for k = 1:length(mpc.user_constraints.nle)
0484 nlc = mpc.user_constraints.nle{k};
0485 fcn = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
0486 hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
0487 om.add_nln_constraint(nlc{1:2}, 1, fcn, hess, nlc{5});
0488 end
0489 end
0490 if nnli
0491 for k = 1:length(mpc.user_constraints.nli)
0492 nlc = mpc.user_constraints.nli{k};
0493 fcn = eval(['@(x)' nlc{3} '(x, nlc{6}{:})']);
0494 hess = eval(['@(x, lam)' nlc{4} '(x, lam, nlc{6}{:})']);
0495 om.add_nln_constraint(nlc{1:2}, 0, fcn, hess, nlc{5});
0496 end
0497 end
0498 if nw
0499 user_cost.N = mpc.N;
0500 user_cost.Cw = Cw;
0501 if ~isempty(fparm)
0502 user_cost.dd = fparm(:, 1);
0503 user_cost.rh = fparm(:, 2);
0504 user_cost.kk = fparm(:, 3);
0505 user_cost.mm = fparm(:, 4);
0506 end
0507 if ~isempty(H)
0508 user_cost.H = H;
0509 end
0510 om.add_legacy_cost('usr', user_cost, user_vars);
0511 end
0512
0513
0514 om = run_userfcn(userfcn, 'formulation', om, mpopt);
0515
0516
0517 cp = om.params_legacy_cost();
0518 [N, H, Cw, rh, mm] = deal(cp.N, cp.H, cp.Cw, cp.rh, cp.mm);
0519 [nw, nx] = size(N);
0520 if nw
0521 if any(cp.dd ~= 1) || any(cp.kk)
0522 if dc
0523 if any(cp.dd ~= 1)
0524 error('opf_setup: DC OPF can only handle legacy user-defined costs with d = 1');
0525 end
0526 if any(cp.kk)
0527 error('opf_setup: DC OPF can only handle legacy user-defined costs with no "dead zone", i.e. k = 0');
0528 end
0529 elseif ~legacy_formulation
0530
0531 legacy_cost_fcn = @(x)opf_legacy_user_cost_fcn(x, cp);
0532 om.add_nln_cost('usr', 1, legacy_cost_fcn);
0533 end
0534 else
0535
0536 if dc || ~legacy_formulation
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547
0548
0549 M = sparse(1:nw, 1:nw, mm, nw, nw);
0550 MN = M * N;
0551 MR = M * rh;
0552 HMR = H * MR;
0553 HtMR = H' * MR;
0554 Q = MN' * H * MN;
0555 c = full(MN' * (Cw - 1/2*(HMR+HtMR)));
0556 k = (1/2 * HtMR - Cw)' * MR;
0557 om.add_quad_cost('usr', Q, c, k);
0558 end
0559 end
0560 end