fm_build.m

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

function fm_build
%FM_BUILD build new component functions (Symbolic Toolbox is needed)
%
%FM_BUILD
%
%see also FM_MAKE FM_COMPONENT
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    19-Dec-2003
%Version:   1.0.1
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Comp Settings Fig Path
global Algeb Buses Initl Param Servc State

% ********************************************************************************
% some control variables
error_v = [];
lasterr('');
null = '0';

% useful strings
c_name = Comp.name;
c_name(1) = upper(c_name(1));

% variable arrays
state_vect = varvect(State.name, ' ');
algeb_vect = varvect(Algeb.name, ' ');
param_vect = varvect(Param.name, ' ');
initl_vect = varvect(Initl.name, ' ');
servc_vect = varvect(Servc.name, ' ');
pq_Servc = 0;

% equation arrays
servc_eq = varvect(Servc.eq,'; ');
state_eq = varvect(State.eq,'; ');
algeb_eq = varvect(Algeb.eq,'; ');

% ********************************************************************************
% check equations
if State.neq > 0
  state_check = strmatch('null',State.eq,'exact');
  if state_check
    error_v = [error_v; ...
               strcat('Differential equation for "', ...
                      State.eqidx(state_check), ...
                      '" has not been defined.')];
  end
end
if Servc.neq > 0
  servc_foo = strcat(Servc.type,Servc.eq);
  servc_check = strmatch('Innernull',servc_foo,'exact');
  servc_check = [servc_check; strmatch('Outputnull',servc_foo,'exact')];
  if servc_check
    error_v = [error_v; ...
               strcat('Service equation for "', ...
                      Servc.eqidx(servc_check), ...
                      '" has not been defined.')];
  end
end

% ********************************************************************************
% check variable usage
servc_idx = [strmatch('Inner',Servc.type); strmatch('Output',Servc.type)];
total_var = [State.name; Algeb.name; Servc.eqidx(servc_idx); Param.name; Initl.name];
total_eqn = [' ',servc_eq,' ',state_eq,' ',algeb_eq,' ',varvect(State.time,'*'),' '];
for i = 1:length(total_var)
  idx = findstr(total_eqn,total_var{i});
  if isempty(idx)
    error_v{end+1,1} = ['The variable "',total_var{i},'" is not used in any equation.'];
  else
    before = total_eqn(idx-1);
    after  = total_eqn(idx+length(total_var{i}));
    check = 1;
    for j = 1:length(idx)
      a = double(after(j));       b = double(before(j));
      a1 = ~isletter(after(j));   a2 = (a ~= 95);  a3 = (a > 57 | a < 48);
      b1 = ~isletter(before(j));  b2 = (b ~= 95);  b3 = (b > 57 | b < 48);
      if a1 & a2 & a3 & b1 & b2 & b3, check = 0; break, end
    end
    if check
      error_v{end+1,1} = ['The variable "',total_var{i}, ...
                          '" is not used in any equation.'];
    end
  end
end

% ********************************************************************************
% symbolic variables
try
  if state_vect, eval(['syms ',state_vect]), end
  if algeb_vect, eval(['syms ',algeb_vect]), end
  if param_vect, eval(['syms ',param_vect]), end
  if servc_vect, eval(['syms ',servc_vect]), end
  if initl_vect, eval(['syms ',initl_vect]), end

  % compute Jacobians matrices (Maple Symbolic Toolbox)
  if ~isempty(state_eq)
    if ~isempty(state_vect)
      eval(['Fx = jacobian([',state_eq, '],[', state_vect, ']);']);
    end
    if ~isempty(algeb_vect)
      eval(['Fy = jacobian([',state_eq, '],[', algeb_vect, ']);']);
    end
    if ~isempty(servc_vect)
      eval(['Fz = jacobian([',state_eq, '],[', servc_vect, ']);']);
    end
    if pq_Servc
      eval(['Fpq = jacobian([',state_eq, '],[', pq_vect, ']);']);
    end
  end
  if ~isempty(algeb_eq)
    if ~isempty(state_vect)
      eval(['Gx = jacobian([',algeb_eq, '],[', state_vect, ']);']);
    end
    if ~isempty(algeb_vect)
      eval(['Gy = jacobian([',algeb_eq, '],[', algeb_vect, ']);']);
    end
    if ~isempty(servc_vect)
      eval(['Gz = jacobian([',algeb_eq, '],[', servc_vect, ']);']);
    end
    if pq_Servc
      eval(['Gpq = jacobian([',alg_eqeb, '],[', pq_vect, ']);']);
    end
  end
  if ~isempty(servc_eq)
    if ~isempty(state_vect)
      eval(['Zx = jacobian([',servc_eq, '],[', state_vect, ']);']);
    end
    if ~isempty(algeb_vect)
      eval(['Zy = jacobian([',servc_eq, '],[', algeb_vect, ']);']);
    end
    if ~isempty(servc_vect)
      eval(['Zz = jacobian([',servc_eq, '],[', servc_vect, ']);']);
    end
    if pq_Servc
      eval(['Zpq = jacobian([',servc_eq, '],[', pq_vect, ']);']);
    end
  end
end

% ********************************************************************************
% check synthax of equations
for i = 1:State.neq
  try
    eval([State.eq{i,1},';'])
  catch
    error_v{end+1,1} = [lasterr, ' (In differential equation "', ...
                        State.eq{i,1}, '")'];
  end
  try
    eval([State.init{i,1},';'])
  catch
    error_v{end+1,1} = [lasterr, ' (In state variable  "', ...
                        State.name{i,1}, '" initialization expression)'];
  end
  state_init{i,1} = vectorize(State.init{i,1});
end
for i = 1:Algeb.neq
  try
    eval([Algeb.eq{i,1},';'])
  catch
    error_v{end+1,1} = [lasterr, ' (In algebraic equation "', ...
                        State.eq{i,1}, '")'];
  end
end
for j = 1:Servc.neq
  try
    eval([Servc.eq{i,1},';'])
  catch
    error_v{end+1,1} = [lasterr, ' (In service equation "', ...
                        Servc.eqidx{i,1}, '")'];
  end
end

% check component name
if isempty(Comp.name)
  error_v{end+1,1} = 'Component name is empty.';
  return
end

% ********************************************************************************
% display errors
if ~isempty(error_v)
  error_v = strcat('Error#',num2str([1:length(error_v)]'),': ',error_v);
  error_v = [{['REPORT OF ERRORS ENCOUNTERED WHILE BUILDING ', ...
               'NEW COMPONENT "',Comp.name,'.m"']}; error_v];
  error_v{end+1,1} = ['BUILDING NEW COMPONENT FILE "', ...
                      Comp.name,'.m" FAILED'];
  fm_disp
  fm_disp(error_v{1:end-1})
  fm_disp(['BUILDING NEW COMPONENT FILE "',Comp.name,'.m" FAILED'])

  fm_update
  set(findobj(Fig.update,'Tag','Listbox1'),'String',error_v, ...
                    'BackgroundColor','w', ...
                    'ForegroundColor','r', ...
                    'Enable','inactive', ...
                    'max',2, ...
                    'Value',[]);
  set(findobj(Fig.update,'Tag','Pushbutton2'),'Enable','off');
  return
end

% ********************************************************************************
% check for previous versions
a = what(Path.psat);
olderfile = strmatch(['fm_',Comp.name,'.m'],a.m,'exact');
if ~isempty(olderfile)
  uiwait(fm_choice(['Overwrite Existing File "fm_',Comp.name,'.m" ?']));
  if ~Settings.ok, return, end
end

% ********************************************************************************
% open new component file
fid = fopen([Path.psat, 'fm_', Comp.name,'.m'], 'wt');
if fid == -1
  fm_disp(['Cannot open file fm_',Comp.name,'. Check permissions'])
  return
end
fprintf(fid, ['function  fm_', Comp.name, '(flag)']);

% write help of the function
if isempty(Comp.descr)
  Comp.descr = ['Algebraic Differential Equation ', ...
                Comp.name, '.m'];
end
fprintf(fid, ['\n\n%%FM_', upper(Comp.name),' defines ',Comp.descr]);

% ********************************************************************************
% data format .con
fprintf(fid, ['\n%%\n%%Data Format ', c_name, '.con:']);
fprintf(fid, '\n%%  col #%d: Bus %d number',[1:Buses.n;1:Buses.n]);
idx_inn = strmatch('Inner',  Servc.type, 'exact');
idx_inp = strmatch('Input',  Servc.type, 'exact');
idx_out = strmatch('Output', Servc.type, 'exact');

fprintf(fid, '\n%%  col #%d: Power rate [MVA]',Buses.n+1);
fprintf(fid, '\n%%  col #%d: Bus %d Voltage Rate [kV]', ...
        [Buses.n+1+[1:Buses.n];1:Buses.n]);
fprintf(fid, '\n%%  col #%d: Frequency rate [Hz]',2*Buses.n+2);
inip = 2*Buses.n+3;
endp  = 2*Buses.n+2+Param.n;
pidx = inip:endp;
for i=1:length(pidx)
  fprintf(fid, '\n%%  col #%d: %s %s [%s]', ...
          pidx(i),Param.name{i},Param.descr{i},Param.unit{i});
end

x_max = [1:State.n];
if State.n
  x_idx = strmatch('None',State.limit(:,1),'exact');
  x_max(x_idx) = [];
end
x_min = [1:State.n];
if State.n
  x_idx = strmatch('None',State.limit(:,2),'exact');
  x_min(x_idx) = [];
end
n_xmax = length(x_max); n_xmin = length(x_min);
for i=1:n_xmax
  fprintf(fid,'\n%%  col #%d: %s',endp+i, ...
          State.limit{x_max(i),1});
end
for i=1:n_xmin
  fprintf(fid,'\n%%  col #%d: %s',endp+n_xmax+i, ...
          State.limit{x_min(i),1});
end

s_max = [1:Servc.neq];
if Servc.n
  s_idx = strmatch('None',Servc.limit(:,1),'exact');
  s_max(s_idx) = [];
end
s_min = [1:Servc.neq];
if Servc.n
  s_idx = strmatch('None',Servc.limit(:,2),'exact');
  s_min(s_idx) = [];
end
n_smax = length(s_max); n_smin = length(s_min);
for i=1:n_smax
  fprintf(fid,'\n%%  col #%d: %s', ...
          endp+n_xmax+n_xmin+i,Servc.limit{s_max(i),1});
end
for i=1:n_smin
  fprintf(fid,'\n%%  col #%d: %s', ...
          endp+n_xmax+n_xmin+n_smax+i,Servc.limit{s_min(i),1});
end

okdata = 0;
nidx = 0;
if ~isempty(idx_inn) | ~isempty(idx_out)
  okdata = 1;
end
if Initl.n | okdata
  fprintf(fid, ['\n%% \n%%Data Structure: ', c_name, '.dat:']);
end
for i=1:Initl.n
  fprintf(fid,'\n%%  col #%d: %s', i,Initl.name{i});
end
if okdata
  nidx = length(idx_inn)+length(idx_out);
  iidx = [idx_inn;idx_out];
  for i=1:nidx
    fprintf(fid,'\n%%  col #%d: %s', ...
            Initl.n+i,Servc.eqidx{iidx(i)});
  end
end

% function calls
fprintf(fid, ['\n%% \n%%FM_', upper(Comp.name),'(FLAG)']);
if Comp.init
  fprintf(fid, ['\n%%   FLAG = 0 -> initialization']);
end
if ~isempty(algeb_eq)
  fprintf(fid, ['\n%%   FLAG = 1 -> algebraic equations']);
  fprintf(fid, ['\n%%   FLAG = 2 -> algebraic Jacobians']);
end
if ~isempty(state_eq);
  fprintf(fid, ['\n%%   FLAG = 3 -> differential equations']);
  fprintf(fid, ['\n%%   FLAG = 4 -> state Jacobians']);
end
if n_xmax | n_xmin > 0
  fprintf(fid, ['\n%%   FLAG = 5 -> non-windup limiters)']);
end
fprintf(fid, '\n%% \n%%Author:    File automatically generated by PSAT');
fprintf(fid, '\n%%Date:      %s',date);

% global variables
fprintf(fid, ['\n\nglobal ',c_name,' DAE Bus Settings']);

% ********************************************************************************
% general settings
fprintf(fid, '\n');
for i=1:State.n
  fprintf(fid,'\n%s = DAE.x(%s.%s);',State.name{i},c_name, ...
          State.name{i});
end
if Algeb.n
  idx_v = strmatch('V',Algeb.name);
  idx_a = strmatch('t',Algeb.name);
  if idx_v
    num_v = strrep(Algeb.name(idx_v),'V','');
    if Buses.n == 1
      fprintf(fid,'\n%s = DAE.V(%s.bus);', ...
              Algeb.name{idx_v},c_name);
    else
      for i=1:length(idx_v)
        fprintf(fid,'\n%s = DAE.V(%s.bus%s);', ...
                Algeb.name{idx_v(i)},c_name,num_v{i});
      end
    end
  end
  if idx_a
    num_a = strrep(Algeb.name(idx_a),'theta','');
    if Buses.n == 1
      fprintf(fid,'\n%s = DAE.a(%s.bus);', ...
              Algeb.name{idx_a},c_name);
    else
      for i=1:length(idx_a)
        fprintf(fid,'\n%s = DAE.a(%s.bus%s);', ...
                Algeb.name{idx_a(i)},c_name,num_a{i});
      end
    end
  end
end
for i=1:Param.n
  fprintf(fid,'\n%s = %s.con(:,%d);',Param.name{i},c_name,pidx(i));
end
for i=1:n_xmax,
  fprintf(fid,'\n%s = %s.con(:,%d);',State.limit{x_max(i),1}, ...
	  c_name,endp+i);
end
for i=1:n_xmin,
  fprintf(fid,'\n%s = %s.con(:,%d);',State.limit{x_min(i),2}, ...
	  c_name,endp+n_xmax+i);
end
for i=1:n_smax,
  fprintf(fid,'\n%s = %s.con(:,%d);',Servc.limit{s_max(i),1}, ...
	  c_name,endp+n_xmax+n_xmin+i);
end
for i=1:n_smin,
  fprintf(fid,'\n%s = %s.con(:,%d);',Servc.limit{s_min(i),2}, ...
	  c_name,endp+n_xmax+n_xmin+n_smax+i);
end
for i=1:Initl.n,
  fprintf(fid,'\n%s = %s.dat(:,%d);',Initl.name{i},c_name,i);
end
for i=1:nidx,
  fprintf(fid,'\n%s = %s.dat(:,%d);',Servc.eqidx{iidx(i)},c_name, ...
	  Initl.n+i);
end

% ********************************************************************************
% initialization
if Comp.init
  fprintf(fid, '\n\nswitch flag\n case 0 %% initialization');
  msg = ['Component'];
  idx_T = [1:State.n];
  idx = strmatch('None',State.time,'exact');
  idx_T(idx) = [];
  if idx_T
    fprintf(fid,'\n\n  %%check time constants');
  end