fmincopf.m

matpower软件下

  A = [ Avl  sparse(size(Avl,1), size(Au,2)-size(Avl,2));
        Au ];
  l = [ lvl;
        lbu ];
  u = [ uvl;
        ubu ];
else               % no z variables
  if (ncony > 0 )  % ... but some y variables from pwl costs; we supply linear
    A = [ Au;      % cost for y variables in last row.
          Avl;
          Ay; 
        sparse(ones(1,ny), ybas:yend, ones(1,ny), 1, yend ) ]; % "linear" cost
    l = [ lbu; 
          lvl;
         -1e10*ones(ncony+1, 1) ];
    u = [ ubu;
          uvl;
          by;
          1e10];
  else                 % No y variables (no pwl costs) and no z variables
    A = [ Au;  Avl ];  % but perhaps we have user linear constraints in Au
    l = [ lbu; lvl ];  % on (theta,V,Pg,Qg) variables.
    u = [ ubu; uvl ];
  end
end

% So, can we do anything good about lambda initialization?
if all(bus(:, LAM_P) == 0)
  bus(:, LAM_P) = (10)*ones(nb, 1);
end

% total number of variables
nxyz = nx+ny+nz;


% --------------------------------------------------------------
% Up to this point, the setup is MINOS-like.  We now adapt
% things for fmincon.

j = sqrt(-1);

% Form a vector with basic info to pass on as a parameter
parms = [ ...
    nb  ;% 1
    ng  ;% 2
    nl  ;% 3
    ny  ;% 4
    nx  ;% 5
    nvl ;% 6
    nz  ;% 7
    nxyz;% 8
    thbas;% 9
    thend;% 10
    vbas;% 11
    vend;% 12
    pgbas;% 13
    pgend;% 14
    qgbas;% 15
    qgend;% 16
    ybas;% 17
    yend;% 18
    zbas;% 19
    zend;% 20
    pmsmbas;% 21
    pmsmend;% 22
    qmsmbas;% 23
    qmsmend;% 24
    sfbas;% 25
    sfend;% 26
    stbas;% 27
    stend;% 28
];

% If there are z variables or y variables the last row of A should be
% holding a linear cost vector of length nxyz.  Let us excise it from A
% explicitly if it exists; otherwise it is zero.
if ny+nz > 0
  nn = size(A,1);
  ccost = full(A(nn, :));
  A(nn, :) = [];
  l(nn) = [];
  u(nn) = [];
else
  ccost = zeros(1, nxyz);
end

% Divide l <= A*x <= u into less than, equal to, greater than, doubly-bounded
% sets.
ieq = find( abs(u-l) <= eps);
igt = find( u >= 1e10);  % unlimited ceiling
ilt = find( l <= -1e10); % unlimited bottom
ibx = find( (abs(u-l) > eps) & (u < 1e10) & (l > -1e10));
Af  = [ A(ilt, :); -A(igt, :); A(ibx, :); -A(ibx, :) ];
bf  = [ u(ilt);   -l(igt);     u(ibx);    -l(ibx)];
Afeq = A(ieq, :);
bfeq = u(ieq);

% bounds on optimization vars; y and z vars unbounded at box bounds;
% if needed, user must do this via the Ax < l mechanism if needed and hope that
% fmincon handles singleton rows elegantly.
UB = Inf * ones(nxyz, 1);
LB = -UB;
LB(thbas+ref-1) = bus(ref, VA)*pi/180;  UB(thbas+ref-1) = bus(ref, VA)*pi/180;
LB(vbas:vend)   = bus(:, VMIN);         UB(vbas:vend)   = bus(:, VMAX);
LB(pgbas:pgend) = gen(:, PMIN)/baseMVA; UB(pgbas:pgend) = gen(:, PMAX)/baseMVA;
LB(qgbas:qgend) = gen(:, QMIN)/baseMVA; UB(qgbas:qgend) = gen(:, QMAX)/baseMVA;

% Compute initial vector
x0 = zeros(nxyz, 1);
x0(thbas:thend) = bus(:, VA) * pi/180;
x0(vbas:vend)   = bus(:, VM);
x0(vbas+gen(:,GEN_BUS)-1) = gen(:, VG);   % buses w. gens init V from gen data
x0(pgbas:pgend) = gen(:, PG) / baseMVA;
x0(qgbas:qgend) = gen(:, QG) / baseMVA;
% no ideas to initialize z, y variables, though, and no mechanism yet
% to ask for user-provided initial z, y.

% build admittance matrices
[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);

% Tolerances
if mpopt(19) == 0   % # iterations
  mpopt(19) = 150 + 2*nb;
end

% basic optimset options needed for fmincon
% fmoptions = optimset('GradObj', 'on', 'Hessian', 'on', 'LargeScale', 'on', ...
%                    'GradConstr', 'on');
fmoptions = optimset('GradObj', 'on', 'LargeScale', 'off', 'GradConstr', 'on');
fmoptions = optimset(fmoptions, 'MaxIter', mpopt(19), 'TolCon', mpopt(16) );
fmoptions = optimset(fmoptions, 'TolX', mpopt(17), 'TolFun', mpopt(18) );
fmoptions.MaxFunEvals = 4 * fmoptions.MaxIter;
if mpopt(31) == 0,
  fmoptions.Display = 'off';
else
  fmoptions.Display = 'iter';
end

Af = full(Af);
Afeq = full(Afeq);
[x, f, info, Output, Lambda, Jac] = ...
  fmincon('costfmin', x0, Af, bf, Afeq, bfeq, LB, UB, 'consfmin', fmoptions, ...
         baseMVA, bus, gen, gencost, branch, areas, Ybus, Yf, Yt, mpopt, ...
         parms, ccost);
success = (info > 0);

% Unpack optimal x
bus(:, VA) = x(thbas:thend)*180/pi;
bus(:, VM) = x(vbas:vend);
gen(:, PG) = baseMVA * x(pgbas:pgend);
gen(:, QG) = baseMVA * x(qgbas:qgend);
gen(:, VG) = bus(gen(:, GEN_BUS), VM);

% reconstruct voltages
Va = x(thbas:thend);
Vm = x(vbas:vend);
V = Vm .* exp(j*Va);

%% compute branch injections
Sf = V(branch(:, F_BUS)) .* conj(Yf * V);  %% cplx pwr at "from" bus, p.u.
St = V(branch(:, T_BUS)) .* conj(Yt * V);  %% cplx pwr at "to" bus, p.u.
branch(:, PF) = real(Sf) * baseMVA;
branch(:, QF) = imag(Sf) * baseMVA;
branch(:, PT) = real(St) * baseMVA;
branch(:, QT) = imag(St) * baseMVA;

% Put in Lagrange multipliers
gen(:, MU_PMAX)  = Lambda.upper(pgbas:pgend) / baseMVA;
gen(:, MU_PMIN)  = Lambda.lower(pgbas:pgend) / baseMVA;
gen(:, MU_QMAX)  = Lambda.upper(qgbas:qgend) / baseMVA;
gen(:, MU_QMIN)  = Lambda.lower(qgbas:qgend) / baseMVA;
bus(:, LAM_P)    = Lambda.eqnonlin(1:nb) / baseMVA;
bus(:, LAM_Q)    = Lambda.eqnonlin(nb+1:2*nb) / baseMVA;
bus(:, MU_VMAX)  = Lambda.upper(vbas:vend);
bus(:, MU_VMIN)  = Lambda.lower(vbas:vend);
branch(:, MU_SF) = Lambda.ineqnonlin(1:nl) / baseMVA; 
branch(:, MU_ST) = Lambda.ineqnonlin(nl+1:2*nl) / baseMVA;

% We are done with standard opf but we may need to provide the
% constraints and their Jacobian also.
if nargout > 7
  [g, geq, dg, dgeq] = consfmin(x, baseMVA, bus, gen, gencost, branch, areas,...
                                Ybus, Yf, Yt, mpopt, parms, ccost);
  g = [ geq; g];
  jac = [ dgeq'; dg']; % true Jacobian organization
end

% Go back to original data.
% reorder generators
gen = gen(inv_gen_ord, :);

% convert to original external bus ordering
[bus, gen, branch, areas] = int2ext(i2e, bus, gen, branch, areas);

% Now create output matrices with all lines, all generators, committed and
% non-committed
genout = genorg;
branchout = branchorg;
genout(comgen, : ) = gen;
branchout(onbranch, :)  = branch;
% And zero out appropriate fields of non-comitted generators and lines
if ~isempty(offgen)
  tmp = zeros(length(offgen), 1);
  genout(offgen, PG) = tmp;
  genout(offgen, QG) = tmp;
  genout(offgen, MU_PMAX) = tmp;
  genout(offgen, MU_PMIN) = tmp;
end
if ~isempty(offbranch)
  tmp = zeros(length(offbranch), 1);
  branchout(offbranch, PF) = tmp;
  branchout(offbranch, QF) = tmp;
  branchout(offbranch, PT) = tmp;
  branchout(offbranch, QT) = tmp;
  branchout(offbranch, MU_SF) = tmp;
  branchout(offbranch, MU_ST) = tmp;
end

et = etime(clock,t1);
if  (nargout == 0) & ( success )
  printpf(baseMVA, bus, genout, branchout, f, info, et, 1, mpopt);
end

if nargout, busout = bus; end


return;