fm_gams.m

电力系统的psat

function fm_gams
% FM_GAMS initialize and call GAMS to solve
%         several kind of Market Clearing Mechanisms
%
% FM_GAMS
%
%GAMS settings are stored in the structure GAMS, with
%the following fields:
%
%      METHOD   1 -> simple auction
%               2 -> linear OPF (DC power flow)
%               3 -> nonlinear OPF (AC power flow)
%               4 -> nonlinear VSC-OPF
%               5 -> maximum loading condition
%               6 -> continuation OPF
%
%      TYPE     1 -> single period auction
%               2 -> multi period auction
%               3 -> VSC single period auction
%               4 -> VSC multi period auction
%
%see also FM_GAMS.GMS, FM_GAMSFIG and
%structures CPF and OPF for futher settings
%
%Author:    Federico Milano
%Date:      29-Jan-2003
%Update:    01-Feb-2003
%Update:    06-Feb-2003
%Version:   1.0.2
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global DAE OPF CPF GAMS Bus File clpsat
global Path Settings Snapshot Varname
global PV PQ SW Line jay Varout
global Supply Demand Rmpl Rmpg Ypdp

if ~autorun('PSAT-GAMS Interface')
  return
end

if DAE.n
  fm_disp(['Dynamic data are not supported within the PSAT-GAMS interface.'],2)
  return  
end

if ~Supply.n,
  fm_disp(['Supply data have to be specified before in order to ', ...
	   'run PSAT-GAMS interface'],2)
  return
end

if ~Demand.n,
  if GAMS.basepg & ~clpsat.init
    Settings.ok = 0;
    uiwait(fm_choice(['It is strongly recommended to exclude ' ...
                      'base generator powers. Do you want to do so?']))
    GAMS.basepg = ~Settings.ok;
  end
  noDem = 1;
  %Demand.bus = PQ.bus;
  %Demand.con = PQ.con(:,[1 2 4 5 4 4 4]);
  %Demand.con = [Demand.con,zeros(PQ.n,8)];
  %Demand.n = PQ.n;
  Demand.bus = 1;
  Demand.con = [1,100,1,zeros(1,12)];
  Demand.n = 1;
else
  noDem = 0;
end

length(Snapshot.V);
if ~GAMS.basepl
  buspl = Snapshot(1).Pl;
  busql = Snapshot(1).Ql;
  Bus.Pl(:) = 0;
  Bus.Ql(:) = 0;
  pqcon = PQ.con(:,[4 5]);
  PQ.con(:,[4 5]) = 0;
end
if ~GAMS.basepg
  ploss = Snapshot(1).Ploss;
  Snapshot(1).Ploss = 0;
  buspg = Snapshot(1).Pg;
  busqg = Snapshot(1).Qg;
  Bus.Pg(:) = 0;
  Bus.Qg(:) = 0;
  pvcon = PV.con(:,4);
  PV.con(:,4) = 0;
end

fm_disp
fm_disp(['---------------------------------------------------------' ...
         '--------'])
fm_disp(' PSAT-GAMS Interface')
fm_disp(['---------------------------------------------------------' ...
         '--------'])
fm_disp

tic
method = GAMS.method;
modelstat = 0;
solvestat = 0;

%if method == 6
%  fm_disp(['Continuation OPF will be available in future PSAT ' ...
%           'versions'],2)
%  return
%end

type = GAMS.type;
omega = GAMS.omega;
if GAMS.method == 6 & GAMS.type ~= 1
  fm_disp(['WARNING: Continuation OPF can be run only with Single' ...
	   ' Period Auctions.'])
  fm_disp('Voltage Stability Constrained OPF will be solved.')
  method = 4;
end
if GAMS.method == 6 & GAMS.flow ~= 1
  fm_disp(['WARNING: Continuation OPF can be run only with Current ' ...
           'Limits.'])
  fm_disp('Current limits in transmission lines will be used.')
  GAMS.flow = 1;
end
if GAMS.type == 3 & GAMS.method ~= 4
  fm_disp(['WARNING: Pareto Set Single Period Auction can be run ' ...
           'only for VSC-OPF.'])
  fm_disp( '         Single Period Auction will be solved.')
  fm_disp
  type = 1;
end
if GAMS.type == 3 & length(GAMS.omega) == 1
  fm_disp(['WARNING: Weighting factor is scalar. ', ...
           'Single Period Auction will be solved.'])
  fm_disp
  type = 1;
end
if GAMS.type == 1 & length(GAMS.omega) > 1
  fm_disp(['WARNING: Weighting factor is a vector. ', ...
           'First omega entry will be used.'])
  fm_disp
  omega = omega(1);
end
if ~rem(GAMS.type,2) & ~Rmpg.n
  type = 1;
  fm_disp(['WARNING: No Ramping data were found. ', ...
           'Single Period Auction will be solved.'])
  fm_disp
end
if GAMS.type == 2 & Rmpg.n & isempty(Ypdp.con)
  type = 4;
  fm_disp(['WARNING: No Power Demand Profile was found. Single ' ...
           'Period Auction with Unit Commitment will be solved.'])
  fm_disp
end

% resetting time vector in case of previous time simulations
if type == 3, Varout.t = []; end

switch method
 case 1, fm_disp(' Simple Auction')
 case 2, fm_disp(' Market Clearing Mechanism')
 case 3, fm_disp(' Standard OPF')
 case 4, fm_disp(' Voltage Stability Constrained OPF')
 case 5, fm_disp(' Maximum Loading Condition')
 case 6, fm_disp(' Continuation OPF')
 case 7, fm_disp(' Congestion Management')
end

switch type
 case 1, fm_disp(' Single-Period Auction')
 case 2, fm_disp(' Multi-Period Auction')
 case 3, fm_disp(' Pareto Set Single-Period Auction')
 case 4, fm_disp(' Single-Period Auction with Unit Commitment')
end

if (GAMS.flatstart | isempty(Snapshot)) & GAMS.method > 2
  DAE.V = ones(Bus.n,1);
  DAE.a = zeros(Bus.n,1);
else
  DAE.V = Snapshot(1).V;
  DAE.a = Snapshot(1).ang;
end

% ------------------------------------------------------------
% Parameter definition
% ------------------------------------------------------------

% dimensions
nBus = int2str(Bus.n);
nPs = int2str(Supply.n);
nPd = int2str(Demand.n);
nQg = int2str(PV.n+SW.n);
nLine = int2str(Line.n);
nBusref = int2str(SW.bus);
nSW = int2str(SW.n);
nPV = int2str(PV.n);
if type == 2
  nH = ['H',num2str(Ypdp.len)];
elseif type == 4
  nH = 'H1';
end

% indexes
iBPs = Supply.bus;
iBPd = Demand.bus;
iBQg = [SW.bus; PV.bus];
if Settings.octave
  Ps_idx = zeros(Bus.n,Supply.n);
  Pd_idx = zeros(Bus.n,Demand.n);
  SW_idx = zeros(Bus.n,SW.n);
  PV_idx = zeros(Bus.n,PV.n);
else
  Ps_idx = sparse(Bus.n,Supply.n);
  Pd_idx = sparse(Bus.n,Demand.n);
  SW_idx = sparse(Bus.n,SW.n);
  PV_idx = sparse(Bus.n,PV.n);
end

for i = 1:Supply.n,
  Ps_idx(Supply.bus(i),i)=1;
end
for i = 1:Demand.n,
  Pd_idx(Demand.bus(i),i)=1;
end
for i = 1:SW.n
  SW_idx(SW.bus(i),i)=1;
end
for i = 1:PV.n
  PV_idx(PV.bus(i),i)=1;
end

% Flows on transmission lines
tps = Line.con(:,11).*exp(jay*Line.con(:,12)*pi/180);
VV = DAE.V.*exp(jay*DAE.a);
r = Line.con(:,8);
rx = Line.con(:,9);
chrg = Line.con(:,10)/2;
z = r + jay*rx;
y = 1./z;
g = real(y);
b = imag(y);
nl = [1:Line.n];
Fij = VV(Line.from).*conj(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
      VV(Line.from).*(jay*chrg))./(tps.*conj(tps)));
Pij0 = real(Fij);
Qij0 = imag(Fij);
Fji = VV(Line.to).*conj((VV(Line.to) - VV(Line.from)./tps).*y + ...
      VV(Line.to).*(jay*chrg));
Pji0 = real(Fji);
Qji0 = imag(Fji);

% Simplified Jacobian Matrices of flows on transmission lines
m = tps.*conj(tps);

y1 = tps.*y./m;
g1 = real(y1);
b1 = imag(y1);
g0 = g./m;
if GAMS.flow == 1 | GAMS.flow == 3
  b0 = chrg;
else
  b0 = (chrg-b)./m;
end
if method == 2
  Li = sparse(1:Line.n,Line.from,b,Line.n,Bus.n);
  Lj = sparse(1:Line.n,Line.to,b,Line.n,Bus.n);
else
  Li = sparse(1:Line.n,Line.from,1,Line.n,Bus.n);
  Lj = sparse(1:Line.n,Line.to,1,Line.n,Bus.n);
end
if Settings.octave
  Li = full(Li);
  Lj = full(Lj);
end
Gh = real(Line.Y);
Bh = imag(Line.Y);
if GAMS.method == 4 | GAMS.method == 6 | GAMS.method == 7
  if GAMS.line
    line_orig = Line.con;
    Y_orig = Line.Y;
    Line.con(GAMS.line,[8 9 10]) = [0, 1e40, 0];
    line_cont = Line.con;
    fm_y;
    Ghc = real(Line.Y);
    Bhc = imag(Line.Y);
    Line.Y = Y_orig;
    Line.con = line_orig;
  else
    Ghc = Gh;
    Bhc = Bh;
  end
end

if GAMS.flow
  Pijmax = Line.con(:,12+GAMS.flow);
  idx_no_I_max = find(Pijmax == 0);
  Pijmax(idx_no_I_max) = 999;
  Pjimax = Pijmax;
else
  Pijmax = 999*Settings.mva*ones(Line.n,1);
  Pjimax = 999*Settings.mva*ones(Line.n,1);
end

% Load power factor
tgphi = zeros(Demand.n,1);
idx = find(Demand.con(:,3) ~= 0);
if ~isempty(idx)
  tgphi(idx) = Demand.con(idx,4)./Demand.con(idx,3);
end

% Bid blocks
Csa = Supply.con(:,7)/Settings.mva;
Csb = Supply.con(:,8);
Csc = Settings.mva*Supply.con(:,9);

Dsa = Supply.con(:,10)/Settings.mva;
Dsb = Supply.con(:,11);
Dsc = Settings.mva*Supply.con(:,12);

Cda = Demand.con(:,8)/Settings.mva;
Cdb = Demand.con(:,9);
Cdc = Settings.mva*Demand.con(:,10);

Dda = Demand.con(:,11)/Settings.mva;
Ddb = Demand.con(:,12);
Ddc = Settings.mva*Demand.con(:,13);

if GAMS.method == 7
  if length(Supply.con(1,:)) < 15
    Rgup = Supply.con(:,8);
    Rgdw = Supply.con(:,8);
  else
    Rgup = Supply.con(:,15);
    Rgdw = Supply.con(:,16);
  end
  if length(Demand.con(1,:)) < 16
    Rdup = Demand.con(:,9);
    Rddw = Demand.con(:,9);
  else
    Rdup = Demand.con(:,16);
    Rddw = Demand.con(:,17);
  end
end
Csmin = min(Csb);

if GAMS.loaddir
  Ps0 = Supply.con(:,3);
  Pd0 = Demand.con(:,3);
else
  Ps0 = Supply.con(:,6);
  Pd0 = Demand.con(:,7);
end
Psmax = Supply.con(:,4);
Pdmax = Demand.con(:,5);
Psmin = Supply.con(:,5);
Pdmin = Demand.con(:,6);

% distributed slack bus coefficients       %% added by liulin
ksw = zeros(Bus.n,1);
if SW.n, ksw(SW.bus) = SW.con(:,11); end
kpv = zeros(Bus.n,1);
if PV.n, kpv(PV.bus) = PV.con(:,10); end
ksu = zeros(Supply.n,1);
for i = 1:PV.n
  idx = find(Supply.bus == PV.bus(i));
  if idx, ksu(idx) = PV.con(i,10); end
end
for i = 1:SW.n
  idx = find(Supply.bus == SW.bus(i));
  if idx, ksu(idx) = SW.con(i,11); end
end

% Fixed powers
Pg0 = Bus.Pg;
Pl0 = Bus.Pl;
Ql0 = Bus.Ql;

% Generator reactive powers and associated limits
Qg0 = zeros(Bus.n,1);
Qg0(iBQg) = Bus.Qg(iBQg);
Qgmin = zeros(Bus.n,1);
Qgmin(iBQg) = [SW.con(:,7); PV.con(:,7)];
Qgmax = zeros(Bus.n,1);
Qgmax(iBQg) = [SW.con(:,6); PV.con(:,6)];

% Voltage limits
V0 = DAE.V;
t0 = DAE.a;
Vmin = zeros(Bus.n,1);
Vmax = zeros(Bus.n,1);
Vmin(PQ.bus) = PQ.con(:,7);
Vmax(PQ.bus) = PQ.con(:,6);
Vmin(PV.bus) = PV.con(:,9);
Vmax(PV.bus) = PV.con(:,8);
Vmin(SW.bus) = SW.con(:,9);
Vmax(SW.bus) = SW.con(:,8);
Vmin(find(Vmin == 0)) = 0.2;
Vmax(find(Vmax == 0)) = 1.5;

% ------------------------------------------------------------
% Data structures
% ------------------------------------------------------------

if GAMS.method < 7
  S.val = [Ps0,Psmax,Psmin,Csa,Csb,Csc,Dsa,Dsb,Dsc,ksu];
  D.val = [Pd0,Pdmax,Pdmin,tgphi,Cda,Cdb,Cdc,Dda,Ddb,Ddc];
else
  S.val = [Ps0,Psmax,Psmin,Csa,Csb,Csc,Dsa,Dsb,Dsc,ksu,Rgup,Rgdw];
  D.val = [Pd0,Pdmax,Pdmin,tgphi,Cda,Cdb,Cdc,Dda,Ddb,Ddc,Rdup,Rddw];
end
X.val = [V0,t0,Pg0,Qg0,Pl0,Ql0,Qgmax,Qgmin,Vmax,Vmin,ksw,kpv];
L.val = [g1,b1,g0,b0,Pijmax,Pjimax];
lambda.val = [GAMS.lmin(1);GAMS.lmax(1);GAMS.omega(1);GAMS.line];

if type == 2 | type == 4
  S.val = [S.val, zeros(Supply.n,8)];
  S.val(Rmpg.sup,[5:9,11,12]) = Rmpg.con(:,[9,5,6,3,4,7,8]);
  S.val(:,4) = Csb;
  S.val(:,10) = Supply.con(:,13);
  S.labels = {cellstr(num2str([1:Supply.n]')), ...
              {'Ps0','Psmax','Psmin','Cs', ...
               'suc','mut','mdt','rut','rdt','u0','y0','z0'}};
else
  if GAMS.method < 7
    S.labels = {cellstr(num2str([1:Supply.n]')), ...
                {'Ps0','Psmax','Psmin','Csa','Csb', ...
                 'Csc','Dsa','Dsb','Dsc','ksu'}};
  else
    S.labels = {cellstr(num2str([1:Supply.n]')), ...
                {'Ps0','Psmax','Psmin','Csa','Csb', ...
                 'Csc','Dsa','Dsb','Dsc','ksu','RGup','RGdw'}};
  end
end

if GAMS.method < 7
  D.labels = {cellstr(num2str([1:Demand.n]')), ...
              {'Pd0','Pdmax','Pdmin','tgphi','Cda', ...
               'Cdb','Cdc','Dda','Ddb','Ddc'}};
else
  D.labels = {cellstr(num2str([1:Demand.n]')), ...
              {'Pd0','Pdmax','Pdmin','tgphi','Cda', ...
               'Cdb','Cdc','Dda','Ddb','Ddc','RDup','RDdw'}};