fm_thload.m

电力系统的psat

function  fm_thload(flag)
% FM_THLOAD defines Thermostatically Controlled Loads
%
% FM_THLOAD(FLAG)
%       FLAG = 0 initialization
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%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 Thload DAE Bus PQ

T = DAE.x(Thload.T);
G = DAE.x(Thload.G);
V = DAE.V(Thload.bus);
Kp = Thload.con(:,3);
Ki = Thload.con(:,4);
Tc = Thload.con(:,5);
T1 = Thload.con(:,6);
Ta = Thload.con(:,7);
Tref = Thload.con(:,8);
G_max = Thload.con(:,9);
K1 = Thload.con(:,10);
KL = Thload.con(:,11);

% Initialization Routines:

switch flag
 case 0

  Pl = zeros(Thload.n,1);
  % get powers and change PQ loads
  for i = 1:Thload.n
    idx = find(PQ.bus == Thload.bus(i));
    if isempty(idx)
      fm_disp(['Error: No PQ Load found at Therm. Contr. ', ...
	       'Load at Bus ', int2str(Thload.bus(i))],2)
    else
      Pl(i) = Thload.con(i,2)*PQ.P0(idx)/100;
      PQ.con(idx,4) = PQ.con(idx,4) - Pl(i);
      if PQ.con(idx,4) == 0 & PQ.con(idx,5) == 0,
	PQ.con(idx,:) = [];
	PQ.bus(idx) = [];
	PQ.n = PQ.n - 1;
	PQ.P0(idx) = [];
	PQ.Q0(idx) = [];
      end
    end
  end
  %Pl = Thload.con(:,2).*Bus.Pl(Thload.bus)/100;
  DAE.x(Thload.G) = Pl./V./V;
  G = DAE.x(Thload.G);
  DAE.x(Thload.T) = Tref;
  Thload.con(:,10) = (Tref-Ta)./Pl;

  idx = find(T1==0);
  if idx
    Thload.con(idx,8) = 1200;
    thwarn(idx,[' Found T1 = 0. Default value T1 = 1200 s ', ...
		'will be used.'])
  end

  idx = find(KL<1);
  if idx
    Thload.con(idx,11) = 2;
    thwarn(idx,' Found KL < 1. Default value KL = 2 will be used.')
  end

  % fix G_max
  Thload.con(:,9) = 2*G;
  fm_disp('Initialization of thermostatically controlled loads completed.')

 case 1 % algebraic equations

  DAE.gp = DAE.gp + sparse(Thload.bus,1,G.*V.*V,Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J12 = DAE.J12 + sparse(Thload.bus,Thload.bus, ...
			     2*G.*V,Bus.n,Bus.n);

 case 3 % differential equations

  DAE.f(Thload.T) = (Ta-T+K1.*G.*V.^2)./T1;

  % allow no dynamics
  no_dyn_G = find(Tc == 0);
  Tc(no_dyn_G) = 1;
  DAE.f(Thload.G) = (-Kp.*(Ta-T+K1.*G.*V.^2)./T1+Ki.*(Tref-T))./Tc;
  DAE.f(Thload.G(no_dyn_G)) = 0;

  % non-windup limits
  idx = find(G >= G_max & DAE.f(Thload.G) > 0);
  if idx, DAE.f(Thload.G(idx)) = 0; end
  idx = find(G <= 0 & DAE.f(Thload.G) < 0);
  if idx, DAE.f(Thload.G(idx)) = 0; end
  DAE.x(Thload.G) = min(DAE.x(Thload.G),G_max);
  DAE.x(Thload.G) = max(DAE.x(Thload.G),0);

 case 4 % Jacobians of state variables

  DAE.Fx = DAE.Fx + sparse(Thload.T,Thload.T,-1./T1,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Thload.T,Thload.G,K1.*V.^2./T1,DAE.n,DAE.n);

  % allow no dynamics
  no_dyn_G = find(Tc == 0);
  Tc(no_dyn_G) = 1;
  DAE.Fx = DAE.Fx + sparse(Thload.G,Thload.T, ...
			   (Kp./T1-Ki)./Tc,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Thload.G,Thload.G, ...
			   -Kp.*K1.*V.^2./T1./Tc,DAE.n,DAE.n);

  DAE.Fy = DAE.Fy + sparse(Thload.T,Thload.bus+Bus.n, ...
			   2*K1.*G.*V./T1,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Thload.G,Thload.bus+Bus.n, ...
			   -2*Kp.*K1.*G.*V./T1./Tc,DAE.n,2*Bus.n);

  DAE.Gx = DAE.Gx + sparse(Thload.bus+Bus.n,Thload.G,V.*V,2*Bus.n,DAE.n);

 case 5 % non-windup limiters

  idx = find((G >= G_max | G <= 0) &  DAE.f(Thload.G) == 0);
  if idx
    global Settings
    k = Thload.G(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = -eye(length(idx));
    else
      DAE.Ac(k,k) = -speye(length(idx));
    end
  end

end

% -------------------------------------------------------------------
% function for creating warning messages
function thwarn(idx, msg)
global Thload Varname
fm_disp(strcat('Warning: Thermostatically controlled load #', ...
    int2str(idx),' at bus #',Varname.name{Thload.bus(idx)},msg))