%GENFORM: Help file describing the generalized OPF formulation used by
%         the fmincopf and MINOPF solvers.
%
%*-*-*-*-*
%CONTENTS
%*-*-*-*-*
%
%1: General OPF problem formulation
%2: Problem data transformations & general linear restrictions
%3: Piecewise linear convex cost formulation using constrained cost variables
%4: Dispatchable loads
%5: Additional information on the generalized formulation structure
%
%
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%1. GENERAL OPF PROBLEM FORMULATION
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
%The generalized OPF formulation used by the fmincon and MINOPF solvers
%can be written as follows:
%
%   min  sum( f1i(Pgi) + f2i(Qgi) ) + c * [x; y; z]
%  x,y,z
%
%subject to
%  g(x) <=> 0          (nonlinear constraints: bus power balance equations
%                                              & branch flow limits)
%  l <= A*[x;y;z] <= u (general linear constraints)
%  xmin <= x <= xmax   (variable bounds: voltage limits, generation limits)
%
%where x = [ Theta; V; Pg; Qg ], and (Pg, Qg) are the ng active and ng
%reactive injections into the network and (V,Theta) are the voltage
%magnitudes and voltage phase angles at each of the nb buses in the
%network. The objective is a generator-wise additive function f(Pg,Qg)
%plus a linear term, including (y, z) which are described below. The
%nonlinear constraints g() are grouped as
%
%  g(x) = [ gp(x);    % nb bus real power mismatch constraints
%           gq(x);    % nb bus reactive power mismatch constraints
%           gsf(x);   % nl branch MVA injection upper limits at "from" end
%           gst(x) ]; % nl branch MVA injection upper limits at "to" end
%
%The additional variables y and z can optionally augment the problem to
%allow casting more general constraints and is also the mechanism though
%which piecewise linear convex costs are modeled. The linear cost vector
%c allows the addition of a general linear cost on the whole set [x;y;z]
%of variables. With this setup, mixed polynomial/piecewise linear costs
%can be dealt with, as well as dispatchable or price-sensitive loads
%requiring a constant power factor.
%
%
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%2. PROBLEM DATA TRANSFORMATIONS AND GENERAL LINEAR RESTRICTIONS
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
%To add general linear constraints of ones own, it is necessary to
%understand the standard transformations performed on the input data
%(bus, gen, branch, areas and gencost tables) before the problem is
%solved in order to know where the optimization variables end up in the x
%vector. All of these transformations are reversed after solving the
%problem so that output data is in the right place in the tables.
%
%The first step filters out inactive generators and branches; original
%tables are saved for data output.
%
%  comgen   = find(gen(:,GEN_STATUS) > 0);     % find online generators
%  onbranch = find(branch(:,BR_STATUS) ~= 0);  % find online branches
%  gen      = gen(comgen, :);
%  branch   = branch(onbranch, :);
%
%The second step is a renumbering of the bus numbers in the bus table so
%that the resulting table contains consecutively-numbered buses starting
%from 1:
%
%  [i2e, bus, gen, branch, areas] = ext2int(bus, gen, branch, areas);
%
%where i2e is saved for inverse reordering at the end. Finally,
%generators are further reordered by bus number:
%
%  ng = size(gen,1);                 % number of generators or injections
%  [tmp, igen] = sort(gen(:, GEN_BUS));
%  [tmp, inv_gen_ord] = sort(igen);  % save for inverse reordering at the end
%  gen  = gen(igen, :);
%  if ng == size(gencost,1)          % This is because gencost might have
%    gencost = gencost(igen, :);     % twice as many rows as gen if there
%  else                              % are reactive injection costs.
%    gencost = gencost( [igen; igen+ng], :);
%  end
%
%Having done this, the variables inside the x vector now have the same
%ordering as in the bus, gen tables:
%
%  x = [ Theta ;    % nb bus voltage angles
%          V   ;    % nb bus voltage magnitudes
%          Pg  ;    % ng active power injections (p.u.) (ascending bus order)
%          Qg ];    % ng reactive power injections (p.u.)(ascending bus order)
%
%and the nonlinear constraints have the same order as in the bus, branch
%tables
%
%  g = [ gp;        % nb real power flow mismatches (p.u.)
%        gq;        % nb reactive power flow mismatches (p.u.)
%        gsf;       % nl "from" end apparent power injection limits (p.u.)
%        gst ];     % nl "to" end apparent power injection limits (p.u.)
%
%With this setup, box bounds on the variables are applied as follows: the
%reference angle is bounded above and below with the value specified for
%it in the original bus table. The V section of x is bounded above and
%below with the corresponding values for VMAX and VMIN in the bus table.
%The Pg and Qg sections of x are bounded above and below with the
%corresponding values for PMAX, PMIN, QMAX and QMIN in the gen table. The
%nonlinear constraints are similarly setup so that gp and gq are equality
%constraints (zero RHS) and the limits for gsf, gst are taken from the
%RATE_A column in the branch table.
%
%The following example illustrates how an additional general linear
%constraint can be added to the problem formulation. In the standard
%solution to case9.m, the voltage angle for bus 7 lags the voltage angle
%in bus 2 by 6.09 degrees. Suppose we want to limit that lag to 5 degrees
%at the most. A linear restriction of the form
%
%  Theta(2) - Theta(7) <= 5 degrees
%
%would do the trick. We have nb = 9 buses, ng = 3 generators and nl = 9
%branches. Therefore the first 9 elements of x are bus voltage angles,
%elements 10:18 are bus voltage magnitudes, elements 19:21 are active
%injections corresponding to the generators in buses 1, 2 and 3 (in that
%order) and elements 22:24 are the corresponding reactive injections.
%Note that in this case the generators in the original data already
%appear in ascending bus order, so no permutation with respect to the
%original data is necessary. Going back to the angle restriction, we see
%that it can be cast as
%
%  [ 0 1 0 0 0 0 -1 0 0 zeros(1,nb+ng+ng) ] * x  <= 5 degrees
%
%We can set up the problem as follows:
%
%  A = sparse([1;1], [2;7], [1;-1], 1, 24);
%  l = -Inf;
%  u = 5 * pi/180;
%  mpopt = mpoption('OPF_ALG', 520); % use fmincon w/generalized formulation
%  opf('case9', A, l, u, mpopt)
%
%which indeed restricts the angular separation to 5 degrees. NOTE: In
%this example, the total number of variables is 24, but if there are any
%piecewise linear cost functions, there may be additional "helper"
%variables used by the solver and in that case the number of columns in A
%may need to be larger. The next section describes how this is done. If
%all costs are polynomial, however, no extra variables are needed.
%
%
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%3. PIECEWISE LINEAR CONVEX COST FORMULATION USING CONSTRAINED COST VARIABLES
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
%The generalized formulation allows for an easy way to model any
%piecewise linear costs. Such a cost curve looks like
%
% cost axis
%     ^
%c2   |                           *
%     |                        *
%     |                     *
%c1   |                  *
%     |          *
%c0   |   *
%     |------------------------------->   Pg axis
%         x0            x1         x2
%        Pmin                     Pmax
%
%This non-differentiable cost can be modeled using one helper cost
%variable for each such cost curve and additional restrictions on this
%variable and Pg, one restriction for each segment of the curve. The
%restrictions build a convex "basin" and they are equivalent to saying
%that the cost variable must lie in the epigraph of the cost curve. When
%the cost is minimized, the cost variable will be pushed against this
%basin. If the helper cost variable is y, then the contribution of the
%generators' cost to the total cost is exactly y, and in the above case
%the two restrictions needed are
%
%  1)  y >= m1*(Pg - x0) + c0   (y must lie above the first segment)
%
%  2)  y >= m2*(Pg - x1) + c1   (y must lie above the second segment)
%
%where, m1 and m2 are the slopes of the two segments. Also needed, of
%course, are the box restrictions on Pg:  Pmin <= Pg <= Pmax. The
%additive part of the cost contributed by this generator is y.
%
%In the generalized OPF formulation, the capability to accept general
%linear constraints is used to introduce new y variables (one for each
%piecewise linear cost in the problem) and constraints (one for each cost
%segment in the problem). The function that builds the coefficient matrix
%for the restrictions is makeAy. Because a linear cost on the y variables
%is also required, the last row of the matrix that is passed to the
%solver is expected to contain not some linear restriction coefficients
%but a linear cost vector on [x; y]. In normal use this is done
%automatically inside fmincopf (or mopf when using MINOPF) and the user
%need not worry about this. To incorporate additional linear constraints,
%however, it is necessary to know in advance how many y variables there
%are so that the coefficient matrix for the user's constraints have a
%matching number of columns to multiply [x; y]. The number of y variables
%is equal to the number of piecewise linear cost curves (both active and
%reactive) there are for the generators that are online.
%
%
%*-*-*-*-*-*-*-*-*-*-*-*
%4.0 DISPATCHABLE LOADS
%*-*-*-*-*-*-*-*-*-*-*-*
%
%In general, dispatchable or price-sensitive loads can be modeled as negative
%real power injections with associated costs. The current test is that if
%PMIN < PMAX == 0 for a generator, then it is really a dispatchable
%load. If a load has a demand curve like the following
%
% quantity
%    ^
%    |
%    |
%P2  |-------------+
%    |             |
%    |             |
%P1  |_____________|_____________+
%    |                           |
%    |                           |
%    |---------------------------+---->   price
%                price1        price2
%
%so that it will consume zero if the price is higher than price2, P1 if
%the price is less than price2 but higher than price1, and P2 if the
%price is equal or lower than price1. Considered as a negative injection,
%the desired dispatch is zero if the price is greater than price2, -P1
%if the price is higher than price1 but lower than price2, and -P2 if the
%price is equal to or lower than price1. This suggests the following
%piecewise linear cost curve:
%
%                                      cost
%                                       ^
%         -P2              -P1         0|
%---------------------------------------*-->  power
%                                     * |0
%                      slope=price2 *   |
%                                 *     |
%                               *       |
%                             *         |
%                           *           |
%                      *                |
%                 *                     |
%            * slope=price1             |
%       *                               |
%
%
%Note that this assumes that the demand blocks can be partially
%dispatched or "split"; if the price trigger is reached half-way through
%the block, the load must accept the partial block. Otherwise, accepting
%or rejecting whole blocks really poses a mixed-integer problem, which is
%outside the scope of the current MATPOWER implementation.
%
%When there are dispatchable loads, the issue of reactive dispatch
%arises. If the QMIN/QMAX generation limits for the "negative generator"
%in question are not set to zero, then the algorithm will dispatch the
%reactive injection to the most convenient value. Since this is not
%normal load behavior, in the generalized formulation it is assumed that
%dispatchable loads maintain a constant power factor. The mechanism for
%posing additional general linear constraints is employed to
%automatically include restrictions for these injections to keep the
%ratio of Pg and Qg constant. This ratio is inferred from the values
%of PMIN and either QMIN (for inductive loads) or QMAX (for capacitive
%loads) in the gen table. It is important to set these appropriately,
%keeping in mind that PG is negative and that for normal inductive
%loads QG should also be negative (a positive reactive load is a
%negative reactive injection). The initial values of the PG and QG
%columns of the gen matrix must be consistent with the ratio defined
%by PMIN and the appropriate Q limit.
%
%
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%5. ADDITIONAL INFORMATION ON THE GENERALIZED FORMULATION STRUCTURE 
%*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
%The advanced user may want to incorporate additional linear constraints
%and/or linear costs to the problem. Sometimes new state variables may
%need to be defined (in addition to x and y discussed above). And,
%sometimes the user may also want to include additional linear costs on
%some or all of the variables, but in most cases fmincopf (or mopf) hide
%these steps from the user and create certain costs and restrictions
%automatically. This creates a conflict between ease of use and
%generality of the software. In this section we explain the behavior
%adopted by fmincopf and mopf when the following factors interact:
%
%a) existence of y variables because of piecewise linear costs
%b) user-provided linear restrictions
%   b1) on whatever the existing variables [x; y] are;
%   b2) with A having more columns than there are elements in [x; y],
%       thus creating extra variables z and an overall optimization vector
%       [x; y; z].
%
%If the user does not provide any additional linear restrictions via the
%A, l, u parameters, then internally there are only x and y-type
%variables; if there are y-type variables for modeling piecewise linear
%costs, then some additional constraints will be constructed to model the
%cost segments. If there are dispatchable loads (with either
%polynomial or piecewise linear costs) then additional constraints will
%be included to maintain a constant power factor.
%
%If a user does provide A, l, u parameters to add general linear
%constraints, then
%
%a) if the number of columns in A is the same as the number of variables
%   in [x; y], then any linear cost vector needed for modeling of the
%   piecewise linear costs is constructed automatically and the user is
%   not expected to provide a cost row in A.
%
%b) if the number of columns in A is greater than the number of elements
%   in [x; y], then the user is responsible to provide
%
%   b1) the appropriate linear restrictions on y and the Pg, Qg sections
%       of the x vector to model each of the segments of the piecewise
%       linear costs as described in previously; this includes both the
%       coefficient matrix and the left and right hand sides of the
%       restrictions. The function makeAy can be used for this, but keep
%       in mind that it will return a coefficient matrix that has only as
%       many columns as elements in [x; y]. It must be padded with enough
%       sparse zero columns on the right to make it conformable with
%       [x; y; z] if there are any z variables.
%
%   b2) the appropriate linear cost in the last row of A and this
%       includes ANY NECESSARY COST COEFFICIENTS ON THE y SECTION OF THE
%       COST VECTOR. Note that l and u must still have the same number of
%       elements as A has rows but the last elements in l and u are
%       meaningless; set them to (-large, +large).
%
%Note that this implies that in order to add linear costs on just the
%[x; y] variables, at least one dummy z variable must be created for the
%interface to provide the user with the opportunity to specify the linear
%cost vector.
%
%Finally, to complete this section we include the structure of
%the Jacobian used internally by MINOS; note that the sparse matrix
%returned by the solvers only includes the rows up to the last of the
%nonlinear constraints (those related to the injection limit on the
%"to" side of the branches).
%
%        nb       nb       ng       ng               ny + nz
%    +------------------------------------+-------------------------+
%    |   dgP      dgP      dgP            |                         |
%nb  |   ---      ---      ---       0    |                         |
%    |   dth      dV       dP             |                         |
%    |                                    |                         |
%    |   dgQ      dgQ               dgQ   |                         |
%nb  |   ---      ---       0       ---   |                         |
%    |   dth      dV                dQ    |                         |
%    |                                    | sparse(2*nb+2*nl,ny+nz) |
%    |   dSf      dSf                     |                         |
%nl  |   ---      ---       0        0    |                         | M = 2*nb
%    |   dth      dV                      |                         |   + 2*nl
%    |                                    |                         |   + mA
%    |   dSt      dSt                     |                         |
%nl  |   ---      ---       0        0    |                         |
%    |   dth      dV                      |                         |
%    +------------------------------------+-------------------------+
%    |                                                              |
%mA  |                                A                             |
%    |                                                              |
%    +--------------------------------------------------------------+
%                         N = 2*nb + 2*ng + ny + nz

%   MATPOWER
%   $Id: genform.m,v 1.7 2005/01/25 14:39:47 ray Exp $
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Autonoma de Manizales
%   Copyright (c) 2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.
