fm_gams.m

这是一个很适合研究和学习用的电力系统仿真软件

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-2006 Federico Milano

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

[u,w] = system('gams');
if u
  fm_disp('GAMS is not properly installed on your system.',2)
  return
end

if ~autorun('PSAT-GAMS Interface',0)
  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(['Exclude (recommended) base generator powers?']))
    GAMS.basepg = ~Settings.ok;
  end
  noDem = 1;
  Demand = add(Demand,'dummy');
else
  noDem = 0;
end

length(Snapshot.V);
if ~GAMS.basepl
  buspl = Snapshot(1).Pl;
  busql = Snapshot(1).Ql;
  Bus.Pl(:) = 0;
  Bus.Ql(:) = 0;
  PQ = pqzero(PQ,'all');
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;
  SW = swzero(SW,'all');
  PV = pvzero(PV,'all');
end

fm_disp
fm_disp('---------------------------------------------------------')
fm_disp(' PSAT-GAMS Interface')
fm_disp('---------------------------------------------------------')
fm_disp

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

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);
nQg = int2str(PV.n+SW.n);
nLine = int2str(Line.n);
nBusref = int2str(SW.bus);

[nSW,SW_idx,ksw] = gams(SW);
[nPV,PV_idx,kpv] = gams(PV);
[nPd,Pd_idx,D] = gams(Demand);
[nPs,Ps_idx,S] = gams(Supply,type);

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];

% 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

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

% 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);
[Qgmax,Qgmin] = fm_qlim(0,0,'all');

% Voltage limits
V0 = DAE.V;
t0 = DAE.a;
[Vmax,Vmin] = fm_vlim(1.5,0.2);

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

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];

X.labels = {cellstr(num2str([1:Bus.n]')), ...
	    {'V0','t0','Pg0','Qg0','Pl0','Ql0', ...
	     'Qgmax','Qgmin','Vmax','Vmin','ksw','kpv'}};
L.labels = {cellstr(num2str([1:Line.n]')), ...
	    {'g','b','g0','b0','Pijmax','Pjimax'}};
lambda.labels = {'lmin','lmax','omega','line'};

X.name = 'X';
L.name = 'N';
lambda.name = 'lambda';

if type == 2
  Ch.val = zeros(Ypdp.len+1,1);
  Ch.val = [0;Ypdp.day(:,Ypdp.d)]'*Ypdp.week(Ypdp.w)*Ypdp.year(Ypdp.y)/1e6;
  Ch.val(1) = Ch.val(2);
  Ch.labels = cell(1,Ypdp.len+1); 
  for idx = 1:Ypdp.len+1
    Ch.labels{idx} = ['H',num2str(idx-1)];
  end
  Ch.name = 'Ch';
elseif type == 4
  Ch.val = [1,1];
  Ch.name = 'Ch';
  Ch.labels = {'H0','H1'};
end

% ------------------------------------------------------------
% Launch GAMS solver
% ------------------------------------------------------------

control = int2str(method);
flow = int2str(GAMS.flow);
currentpath = pwd;
file = 'fm_gams';
if ~rem(type,2)
  file = [file,'2'];
end
if GAMS.libinclude
  file = [file,' ',GAMS.ldir];
end

if clpsat.init | ispc, cd(Path.psat), end

switch control

 % ------------------------------------------------------------------
 case '1' % S I M P L E   A U C T I O N
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction
    [Ps,Pd,MCP,modelstat,solvestat] = psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV, ...
                       nBusref,control,flow,S,D,X);
    [Pij,Pji,Qg] = updatePF(Pd,Ps,iBQg);

  elseif ~rem(type,2) % Single/Multi Period Auction with UC
    [Ps,Pd,MCP,modelstat,solvestat] = psatgams(file,nBus,nLine,nPs,nPd,nSW,nPV, ...
                       nBusref,nH,control,flow,S,D,X,Ch);
    numh = size(MCP,1);
    a = zeros(numh,Bus.n);
    V = zeros(numh,Bus.n);
    Qg = zeros(numh,Bus.n);
    Pij = zeros(numh,Line.n);
    Qij = zeros(numh,Line.n);
    for i = 1:numh
      [Piji,Pjii,Qgi] = updatePF(Pd(i,:)',Ps(i,:)',iBQg);
      a(i,:) = DAE.a';
      V(i,:) = DAE.V';
      Pij(i,:) = Piji';
      Pji(i,:) = Pjii';
      Qg(i,iBQg) = Qgi';
    end
    ro = MCP*ones(1,Bus.n);

  end

 % ------------------------------------------------------------------
 case '2' % M A R K E T   C L E A R I N G   M E C H A N I S M
 % ------------------------------------------------------------------

  if type == 1 % Single Period Auction
    [Ps,Pd,MCP,modelstat,solvestat] = psatgams(file,nBus, ...
                       nLine,nPs,nPd,nSW,nPV,nBusref,control, ...
                       flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx, ...
                       S,D,X,L);
    [Pij,Pji,Qg] = updatePF(Pd,Ps,iBQg);

  elseif ~rem(type,2) % Single/Multi Period Auction with UC
    [Ps,Pd,MCP,modelstat,solvestat] = psatgams(file,nBus, ...
                       nLine,nPs,nPd,nSW,nPV,nBusref,nH,control, ...
                       flow,Gh,Bh,Li,Lj,Ps_idx,Pd_idx,SW_idx,PV_idx, ...
                       S,D,X,L,Ch);
    numh = size(MCP,1);
    a = zeros(numh,Bus.n);
    V = zeros(numh,Bus.n);
    Qg = zeros(numh,Bus.n);
    Pij = zeros(numh,Line.n);
    Qij = zeros(numh,Line.n);
    for i = 1:numh