LPconstr

------------------------------  deprecated  ------------------------------
   OPF solvers based on LPCONSTR to be removed in a future version.
--------------------------------------------------------------------------
LPCONSTR  Finds solution of NLP problem based on successive LP.
   The key is to set up the problem as follows:
           Min f(xi, xo)
       S.T. g1(xi, xo) =0
            g2(xi, xo) =<0
   where the number of equations in g1 is the same as the number of
   elements in xi.

   [X, LAMBDA, CONVERGED]=LPCONSTR('FUN',x, mpopt, step0 ,VLB,VUB,'GRADFUN', 
   'LPEQUSVR', P1,P2,..) starts at x and finds a constrained minimum to
   the function which is described in FUN (usually an M-file: FUN.M).
   The function 'FUN' should return two arguments: a scalar value of the
   function to be minimized, F, and a matrix of constraints, G:
   [F,G]=FUN(X). F is minimized such that G < zeros(G).

   LPCONSTR allows a vector of optional parameters to be defined. For 
   more information type HELP LPOPTION.

   VLB,VUB define a set of lower and upper bounds on the design variables, X, 
   so that the solution is always in the range VLB <= X <= VUB.

   The function 'GRADFUN' is entered which returns the partial derivatives 
   of the function and the constraints at X:  [gf,GC] = GRADFUN(X).

   The problem-dependent parameters P1,P2,... directly are passed to the 
   functions FUN and GRADFUN: FUN(X,P1,P2,...) and GRADFUN(X,P1,P2,...).

   LAMBDA contains the Lagrange multipliers.

   to be worked out:
   write a generalizer equation solver

   See also LPOPF_SOLVER.

0001 function [x, lambda, converged] = LPconstr(FUN,x,mpopt,step0,VLB,VUB,GRADFUN,LPEQUSVR, ...
0002                     P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15)
0003 %------------------------------  deprecated  ------------------------------
0004 %   OPF solvers based on LPCONSTR to be removed in a future version.
0005 %--------------------------------------------------------------------------
0006 %LPCONSTR  Finds solution of NLP problem based on successive LP.
0007 %   The key is to set up the problem as follows:
0008 %           Min f(xi, xo)
0009 %       S.T. g1(xi, xo) =0
0010 %            g2(xi, xo) =<0
0011 %   where the number of equations in g1 is the same as the number of
0012 %   elements in xi.
0013 %
0014 %   [X, LAMBDA, CONVERGED]=LPCONSTR('FUN',x, mpopt, step0 ,VLB,VUB,'GRADFUN',
0015 %   'LPEQUSVR', P1,P2,..) starts at x and finds a constrained minimum to
0016 %   the function which is described in FUN (usually an M-file: FUN.M).
0017 %   The function 'FUN' should return two arguments: a scalar value of the
0018 %   function to be minimized, F, and a matrix of constraints, G:
0019 %   [F,G]=FUN(X). F is minimized such that G < zeros(G).
0020 %
0021 %   LPCONSTR allows a vector of optional parameters to be defined. For
0022 %   more information type HELP LPOPTION.
0023 %
0024 %   VLB,VUB define a set of lower and upper bounds on the design variables, X,
0025 %   so that the solution is always in the range VLB <= X <= VUB.
0026 %
0027 %   The function 'GRADFUN' is entered which returns the partial derivatives
0028 %   of the function and the constraints at X:  [gf,GC] = GRADFUN(X).
0029 %
0030 %   The problem-dependent parameters P1,P2,... directly are passed to the
0031 %   functions FUN and GRADFUN: FUN(X,P1,P2,...) and GRADFUN(X,P1,P2,...).
0032 %
0033 %   LAMBDA contains the Lagrange multipliers.
0034 %
0035 %   to be worked out:
0036 %   write a generalizer equation solver
0037 %
0038 %   See also LPOPF_SOLVER.
0039 
0040 %   MATPOWER
0041 %   $Id: LPconstr.m,v 1.12 2010/04/26 19:45:25 ray Exp $
0042 %   by Deqiang (David) Gan, PSERC Cornell & Zhejiang University
0043 %   Copyright (c) 1996-2010 by Power System Engineering Research Center (PSERC)
0044 %
0045 %   This file is part of MATPOWER.
0046 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0047 %
0048 %   MATPOWER is free software: you can redistribute it and/or modify
0049 %   it under the terms of the GNU General Public License as published
0050 %   by the Free Software Foundation, either version 3 of the License,
0051 %   or (at your option) any later version.
0052 %
0053 %   MATPOWER is distributed in the hope that it will be useful,
0054 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
0055 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0056 %   GNU General Public License for more details.
0057 %
0058 %   You should have received a copy of the GNU General Public License
0059 %   along with MATPOWER. If not, see <http://www.gnu.org/licenses/>.
0060 %
0061 %   Additional permission under GNU GPL version 3 section 7
0062 %
0063 %   If you modify MATPOWER, or any covered work, to interface with
0064 %   other modules (such as MATLAB code and MEX-files) available in a
0065 %   MATLAB(R) or comparable environment containing parts covered
0066 %   under other licensing terms, the licensors of MATPOWER grant
0067 %   you additional permission to convey the resulting work.
0068 
0069 % ------------------------------ setting up -----------------------------------
0070 
0071 if nargin < 8, error('\ LPconstr needs more arguments ! \   '); end
0072 
0073 
0074 nvars = length(x);
0075 nequ =  mpopt(15);
0076 
0077 % set up the arguments of FUN
0078 if ~any(FUN<48) % Check alphanumeric
0079         etype = 1;
0080         evalstr = [FUN,];
0081         evalstr=[evalstr, '(x'];
0082         for i=1:nargin - 8
0083                 etype = 2;
0084                 evalstr = [evalstr,',P',int2str(i)];
0085         end
0086         evalstr = [evalstr, ')'];
0087 else
0088         etype = 3;
0089         evalstr=[FUN,'; g=g(:);'];
0090 end
0091 
0092 %set up the arguments of GRADFUN
0093 if ~any(GRADFUN<48) % Check alphanumeric
0094         gtype = 1;
0095         evalstr2 = [GRADFUN,'(x'];
0096         for i=1:nargin - 8
0097                 gtype = 2;
0098                 evalstr2 = [evalstr2,',P',int2str(i)];
0099         end
0100         evalstr2 = [evalstr2, ')'];
0101 else
0102         gtype = 3;
0103         evalstr2=[GRADFUN,';'];
0104 end
0105 
0106 %set up the arguments of LPEQUSVR
0107 if ~any(LPEQUSVR<48) % Check alphanumeric
0108         lpeqtype = 1;
0109         evalstr3 = [LPEQUSVR,'(x'];
0110         for i=1:nargin - 8
0111                 lpeqtype = 2;
0112                 evalstr3 = [evalstr3,',P',int2str(i)];
0113         end
0114         evalstr3 = [evalstr3, ')'];
0115 else
0116         lpeqtype = 3;
0117         evalstr3=[LPEQUSVR,';'];
0118 end
0119 
0120 % ----------------------------- the main loop ----------------------------------
0121 verbose = mpopt(31);
0122 itcounter = 0;
0123 runcounter = 1;
0124 
0125 stepsize = step0 * 0.02; % use this small stpesize to detect how close to optimum, so to choose better stepsize
0126 
0127 %stepsize = step0;
0128 %fprintf('\n LPconstr does not adaptively choose starting point \n');
0129 
0130 f_best = 9.9e15;
0131 f_best_run = 9.9e15;
0132 max_slackvar_last = 9.9e15;
0133 converged = 0;
0134 
0135 if verbose
0136     fprintf(' it   obj function   max violation  max slack var    norm grad       norm dx\n');
0137     fprintf('----  -------------  -------------  -------------  -------------  -------------\n');
0138 end
0139 while (converged == 0) && (itcounter < mpopt(22)) && (runcounter < mpopt(23))
0140 
0141     itcounter = itcounter + 1; 
0142     if verbose, fprintf('%3d ', itcounter); end
0143 
0144 
0145     % ----- fix xi temporarily, solve equations g1(xi, xo)=0 to get xo (by Newton method).
0146     if lpeqtype == 1,
0147         [x, success_lf] = feval(LPEQUSVR,x);
0148     elseif lpeqtype == 2
0149         [x, success_lf] = eval(evalstr3);
0150     else
0151         eval(evalstr3);
0152     end
0153 
0154 
0155     if success_lf == 0; 
0156         fprintf('\n      Load flow did not converge. LPconstr restarted with reduced stepsize! '); 
0157         x = xbackup;
0158         stepsize = 0.7*stepsize;
0159     end
0160 
0161 
0162     % ----- compute f, g, df_dx, dg_dx
0163 
0164     if etype == 1,              % compute g(x)
0165             [f, g] = feval(FUN,x);
0166     elseif etype == 2
0167             [f, g] = eval(evalstr);
0168     else
0169             eval(evalstr);
0170     end
0171     if gtype == 1               % compute jacobian matrix
0172             [df_dx, dg_dx] = feval(GRADFUN, x);
0173     elseif gtype == 2
0174             [df_dx, dg_dx] = eval(evalstr2);
0175     else
0176             eval(evalstr2);
0177     end
0178     dg_dx = dg_dx';
0179     max_g = max(g);
0180 
0181     if verbose, fprintf('   %-12.6g   %-12.6g', f, max_g); end
0182 
0183 
0184 
0185     % ----- solve the linearized NP, that is, solve a LP to get dx
0186     a_lp = dg_dx; f_lp = df_dx; rhs_lp = -g; vubdx = stepsize; vlbdx = -stepsize;
0187     if isempty(VUB) ~= 1 || isempty(VLB) ~= 1
0188         error('sorry, at this stage LPconstr can not solve a problem with VLB or VUB ');
0189     end
0190 
0191     % put slack variable into the LP problem such that the LP problem is always feasible
0192 
0193     temp = find( ( g./(abs(g) + ones(length(g), 1) ))  > 0.1*mpopt(16));
0194 
0195 
0196     if isempty(temp) ~= 1
0197             n_slack = length(temp);
0198             if issparse(a_lp)
0199                 a_lp = [a_lp, sparse(temp, 1:n_slack, -1, size(a_lp,1), n_slack)];
0200             else
0201                 a_lp = [a_lp, full(sparse(temp, 1:n_slack, -1, size(a_lp,1), n_slack))];
0202             end
0203             vubdx = [vubdx; g(temp) + 1.0e4*ones(n_slack, 1)];
0204             vlbdx = [vlbdx; zeros(n_slack, 1)];
0205             f_lp = [f_lp; 9.9e6 * max(df_dx) * ones(n_slack, 1)];
0206     end
0207 
0208 
0209     % Ray's heuristics of deleting constraints
0210 
0211     if itcounter ==1
0212             idx_workc = [];
0213             flag_workc = zeros(3 * length(rhs_lp) + 2 * nvars, 1);
0214     else
0215             flag_workc = flag_workc - 1;
0216             flag_workc(idx_bindc) = 20 * ones(size(idx_bindc));
0217 
0218             if itcounter > 20
0219                     idx_workc = find(flag_workc > 0);
0220             end
0221     end
0222 
0223 
0224 
0225     [dx, lambda, idx_workc, idx_bindc] = LPsetup(a_lp, f_lp, rhs_lp, nequ, vlbdx, vubdx, idx_workc, mpopt);
0226 
0227 
0228     if length(dx) == nvars
0229         max_slackvar = 0;
0230     else
0231         max_slackvar = max(dx(nvars+1:length(dx))); if max_slackvar < 1.0e-8, max_slackvar = 0; end;
0232     end
0233 
0234     if verbose, fprintf('   %-12.6g', max_slackvar); end
0235 
0236 
0237     dx = dx(1 : nvars); % stripe off the reduendent slack variables
0238 
0239 
0240     % ----- update x, compute the objective function
0241 
0242     x = x + dx;
0243     xbackup = x;
0244     dL_dx = df_dx + dg_dx' * lambda;    % at optimal point, dL_dx should be zero (from KT condition)
0245     %norm_df = norm(df_dx, inf);
0246     norm_dL = norm(dL_dx, inf);
0247     if abs(f) < 1.0e-10
0248         norm_grad = norm_dL;
0249     else    
0250         norm_grad = norm_dL/abs(f);
0251         %norm_grad = norm_dL/norm_df;  % this is more stringent
0252 
0253     end
0254     norm_dx = norm(dx ./ step0, inf);
0255 
0256     if verbose, fprintf('   %-12.6g   %-12.6g\n', norm_grad, norm_dx); end
0257 
0258     % ----- check stopping conditions
0259 
0260     if (norm_grad < mpopt(20)) && (max_g < mpopt(16)) && (norm_dx < mpopt(21))
0261         converged = 1;  break;
0262     end
0263 
0264 %   if max_slackvar > 1.0e-8 && itcounter > 60, break; end
0265 
0266 
0267     if norm_dx < 0.05 * mpopt(21),      % stepsize is overly small, so what is happening?
0268 
0269         if max_g < mpopt(16) && abs(f_best - f_best_run)/f_best_run < 1.0e-4 
0270 
0271             % The solution is the same as that we got in previous run. So we conclude that
0272             % the stopping conditions are overly stringent, and LPconstr HAS found the solution.
0273             converged = 1;  
0274             break;  
0275 
0276         else
0277             % stepsize is overly small to make good progress, we'd better restart using larger stepsize
0278             f_best_run = f_best;
0279             stepsize = 0.4* step0; 
0280 
0281             if verbose
0282                 fprintf('\n----- restarted with larger stepsize\n');
0283             end
0284 
0285             runcounter = runcounter + 1;
0286         end;
0287     end
0288 
0289 
0290     % ----- adjust stepsize
0291 
0292     if itcounter == 1                           % the 1th iteration is a trial one
0293                                         % whihc sets up starting stepsize
0294         if     norm_grad <       mpopt(20)
0295                         stepsize = 0.015 * step0;   % use extra-small stepsize
0296         elseif norm_grad < 2.0 * mpopt(20)
0297                         stepsize = 0.05 * step0;    % use very small stepsize
0298         elseif norm_grad < 4.0 * mpopt(20)
0299                         stepsize = 0.3  * step0;    % use small stepsize
0300         elseif norm_grad < 6.0 * mpopt(20)
0301                         stepsize = 0.6  * step0;    % use less small stepsize
0302         else
0303                         stepsize =     step0;       % use large stepsize
0304         end
0305     end
0306 
0307     if itcounter > 2
0308         if max_slackvar > max_slackvar_last + 1.0e-10
0309             stepsize = 0.7* stepsize; 
0310         end
0311 
0312         if max_slackvar < 1.0e-7        % the trust region method
0313             actual_df  = f_last - f;
0314             if abs(predict_df) > 1.0e-12
0315                     ratio = actual_df/predict_df;
0316             else
0317                 ratio = -99999;
0318             end
0319 
0320             if ratio < 0.25 || f > f_last * 0.9999
0321                     stepsize = 0.5 * stepsize;
0322             elseif ratio > 0.80
0323                     stepsize = 1.05 * stepsize;
0324             end
0325 
0326             if norm(stepsize ./ step0, inf) > 3.0, stepsize = 3*step0; end;    % ceiling of stepsize
0327         end;
0328     end
0329 
0330     max_slackvar_last = max_slackvar;
0331     f_best = min(f, f_best);
0332     f_last = f;
0333     predict_df = -(df_dx(1:nvars))' * dx(1:nvars);
0334 end
0335 
0336 % ------ recompute f and g
0337 if etype == 1,              % compute g(x)
0338             [f, g] = feval(FUN,x);
0339 elseif etype == 2
0340             [f, g] = eval(evalstr);
0341 else
0342         eval(evalstr);
0343 end
0344 
0345 i = find(g < -mpopt(16));
0346 lambda(i) = zeros(size(i));