lpconstr.m

关于电力系统的计算程序!最新的改进法计算潮流程序！用matlab实现的！ 适用于电力行业的相关从业人员！

function [x, lambda, converged] = LPconstr(FUN,x,idx_xi, mpopt,step0,VLB,VUB,GRADFUN,LPEQUSVR, ...
					P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15)
%
%LPCONSTR Finds the solution of a nonlinear programming problem based on 
%successive linear programm. 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

%   MATPOWER Version 2.0
%   by Deqiang (David) Gan, PSERC Cornell    12/22/97
%   Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/ for more info.

% ------------------------------ setting up -----------------------------------

if nargin < 9, error('\ LPconstr needs more arguments ! \   '); end


nvars = length(x);
nequ =  mpopt(15);

% set up the arguments of FUN
if ~any(FUN<48) % Check alphanumeric
        etype = 1;
        evalstr = [FUN,];
        evalstr=[evalstr, '(x'];
        for i=1:nargin - 9
                etype = 2;
                evalstr = [evalstr,',P',int2str(i)];
        end
        evalstr = [evalstr, ')'];
else
        etype = 3;
        evalstr=[FUN,'; g=g(:);'];
end

%set up the arguments of GRADFUN
if ~any(GRADFUN<48) % Check alphanumeric
        gtype = 1;
        evalstr2 = [GRADFUN,'(x'];
        for i=1:nargin - 9
                gtype = 2;
                evalstr2 = [evalstr2,',P',int2str(i)];
        end
        evalstr2 = [evalstr2, ')'];
else
        gtype = 3;
        evalstr2=[GRADFUN,';'];
end

%set up the arguments of LPEQUSVR
if ~any(LPEQUSVR<48) % Check alphanumeric
        lpeqtype = 1;
        evalstr3 = [LPEQUSVR,'(x'];
        for i=1:nargin - 9
                lpeqtype = 2;
                evalstr3 = [evalstr3,',P',int2str(i)];
        end
        evalstr3 = [evalstr3, ')'];
else
        lpeqtype = 3;
        evalstr3=[LPEQUSVR,';'];
end

% ----------------------------- the main loop ----------------------------------
tequationer = 0;
tcomputefg = 0;
tsetuplp = 0;
tsolvelp = 0;
verbose = mpopt(31);
itcounter = 0;
runcounter = 1;

stepsize = step0 * 0.02; % use this small stpesize to detect how close to optimum, so to choose better stepsize

%stepsize = step0;
%fprintf('\n LPconstr does not adaptively choose starting point \n');

f_best = 9.9e15;
f_best_run = 9.9e15;
max_slackvar_last = 9.9e15;
converged = 0;

if verbose
	fprintf(' it   obj function   max violation  max slack var    norm grad       norm dx\n');
	fprintf('----  -------------  -------------  -------------  -------------  -------------\n');
end
while converged ==0 & itcounter < mpopt(22) & runcounter < mpopt(23)

	itcounter = itcounter + 1; 
	if verbose, fprintf('%3d ', itcounter); end


	% ----- fix xi temporarily, solve equations g1(xi, xo)=0 to get xo (by Newton method).

	if isempty(idx_xi) == 1

	        if lpeqtype == 1,
        	        [x, success_lf] = feval(LPEQUSVR,x);
        	elseif lpeqtype == 2
                	[x, success_lf] = eval(evalstr3);
	        else
        	        eval(evalstr3);
	        end

	else
		temp = ones(nvars, 1); temp(idx_xi) = zeros(length(idx_xi), 1);
		idx_xo = find(temp);

		success_lf = 0; counter_lf = 0;
		while success_lf == 0 & counter_lf < 10

			counter_lf = counter_lf + 1;
			if etype == 1,				% compute g(x)
		        	[f, g] = feval(FUN,x);
			elseif etype == 2
        			[f, g] = eval(evalstr);
			else
			        eval(evalstr);
			end

		        if gtype == 1				% compute jacobian matrix
        		        [df_dx, dg_dx] = feval(GRADFUN, x);
        		elseif gtype == 2
                		[df_dx, dg_dx] = eval(evalstr2);
		        else
				eval(evalstr2);
		        end
			dg_dx = sparse(dg_dx');
		
			dxo = - dg_dx(1:nequ, idx_xo) \ g(1:nequ);
			x(idx_xo) = x(idx_xo) + dxo;
			if norm(dxo, inf) < 1.0e-6; success_lf = 1; end
		end

	end


	if success_lf == 0; 
		fprintf('\n      Load flow did not converge. LPconstr restarted with reduced stepsize! '); 
		x = xbackup;
		stepsize = 0.7*stepsize;
	end


	% ----- compute f, g, df_dx, dg_dx 

	if etype == 1,				% compute g(x)
        	[f, g] = feval(FUN,x);
	elseif etype == 2
        	[f, g] = eval(evalstr);
	else
	        eval(evalstr);
	end

        if gtype == 1				% compute jacobian matrix
                [df_dx, dg_dx] = feval(GRADFUN, x);
        elseif gtype == 2
                [df_dx, dg_dx] = eval(evalstr2);
        else
		eval(evalstr2);
        end
	dg_dx = sparse(dg_dx');


	max_g = max(g);

	if verbose, fprintf('   %-12.6g   %-12.6g', f, max_g); end



	% ----- solve the lineaized NP, that is, solve a LP to get dx


	a_lp = dg_dx; f_lp = df_dx; rhs_lp = -g; vubdx = stepsize; vlbdx = -stepsize;
	if isempty(VUB) ~= 1 | isempty(VLB) ~= 1
		error(' sorry, at this stage LPconstr can not solve a problem with VLB or VUB ');
	end

	% put slack variable into the LP problem such that the LP problem is always feasible

        actual_violation = 0;


        temp = find( ( g./(abs(g) + ones(length(g), 1) ))  > 0.1*mpopt(16));


        if isempty(temp) ~= 1
                n_slack = length(temp);
                a_lp = [a_lp, sparse(temp, 1:n_slack, -1, size(a_lp,1), n_slack)];
                vubdx = [vubdx; g(temp) + 1.0e4*ones(n_slack, 1)];
                vlbdx = [vlbdx; zeros(n_slack, 1)];
                f_lp = [f_lp; 9.9e6 * max(df_dx) * ones(n_slack, 1)];
        end


        % Ray's heuristics of deleting constraints

        if itcounter ==1
                idx_workc = [];
                flag_workc = zeros(3 * length(rhs_lp) + 2 * nvars, 1);
        else
                flag_workc = flag_workc - 1;
                flag_workc(idx_bindc) = 20 * ones(size(idx_bindc));

                if itcounter > 20
                        idx_workc = find(flag_workc > 0);
                end
        end;



	[dx, lambda, idx_workc, idx_bindc] = LPsetup(a_lp, f_lp, rhs_lp, nequ, vlbdx, vubdx, idx_workc, mpopt);


	if length(dx) == nvars
		max_slackvar = 0;
	else
		max_slackvar = max(dx(nvars+1:length(dx))); if max_slackvar < 1.0e-8, max_slackvar = 0; end;
	end

	if verbose, fprintf('   %-12.6g', max_slackvar); end


	dx = dx(1 : nvars); % stripe off the reduendent slack variables


	% ----- update x, compute the objective function

	x = x + dx;
	xbackup = x;
	dL_dx = df_dx + dg_dx' * lambda;	% at optimal point, dL_dx should be zero (from KT condition)
	norm_df = norm(df_dx, inf); norm_dL = norm(dL_dx, inf);
	if abs(f) < 1.0e-10
		norm_grad = norm_dL;
	else	
		norm_grad = norm_dL/abs(f);
		%norm_grad = norm_dL/norm_df;  % this is more stringent

	end
	norm_dx = norm(dx ./ step0, inf);

	if verbose, fprintf('   %-12.6g   %-12.6g\n', norm_grad, norm_dx); end

	% ----- check stopping conditions

	if norm_grad < mpopt(20) 			...
	& max_g      < mpopt(16) 			...
	& norm_dx    < mpopt(21), 
		converged = 1;  break; 
	end;

%	if max_slackvar > 1.0e-8 & itcounter > 60, break; end


	if norm_dx < 0.05 * mpopt(21), 		% stepsize is overly small, so what is happening?

		if max_g < mpopt(16) & abs(f_best - f_best_run)/f_best_run < 1.0e-4 

			% The solution is the same as that we got in previous run. So we conclude that
			% the stopping conditions are overly stringent, and LPconstr HAS found the solution.
			converged = 1;	
			break;	

		else
			% stepsize is overly small to make good progress, we'd better restart using larger stepsize
			f_best_run = f_best;
			stepsize = 0.4* step0; 

			if verbose
				fprintf('\n----- restarted with larger stepsize\n');
			end

			runcounter = runcounter + 1;
		end;
	end


	% ----- adjust stepsize

	if itcounter == 1							% the 1th iteration is a trial one
										% whihc sets up starting stepsize
		if     norm_grad <       mpopt(20)
						stepsize = 0.015 * step0;	% use extra-small stepsize
		elseif norm_grad < 2.0 * mpopt(20)
						stepsize = 0.05 * step0;	% use very small stepsize
		elseif norm_grad < 4.0 * mpopt(20)
						stepsize = 0.3  * step0;	% use small stepsize
		elseif norm_grad < 6.0 * mpopt(20)
						stepsize = 0.6  * step0;	% use less small stepsize
		else
						stepsize =     step0;		% use large stepsize			
		end
	end

	if itcounter > 2
                if max_slackvar > max_slackvar_last + 1.0e-10
			stepsize = 0.7* stepsize; 
                end

		if max_slackvar < 1.0e-7        % the trust region method
		        actual_df  = f_last - f;
			if abs(predict_df) > 1.0e-12
		        	ratio = actual_df/predict_df;
			else
				ratio = -99999;
			end

	                if ratio < 0.25 | f > f_last * 0.9999
        	                stepsize = 0.5 * stepsize;
           		elseif ratio > 0.80
        	                stepsize = 1.05 * stepsize;
			end

              		if norm(stepsize ./ step0, inf) > 3.0, stepsize = 3*step0; end;    % ceiling of stepsize
		end;
	end

	max_slackvar_last = max_slackvar;
	f_best = min(f, f_best);
        f_last = f;
        predict_df = -(df_dx(1:nvars))' * dx(1:nvars);
end