fm_build.m

来自「这是一个很适合研究和学习用的电力系统仿真软件」· M 代码 · 共 818 行 · 第 1/2 页

  for i=1:length(idx_T),
    fprintf(fid,['\n  idx = find(%s == 0);\n  if idx\n    ', ...
                 Comp.name,'warn(idx, ''Time constant %s ', ...
                 'cannot be zero. %s = 0.001 s will be used.''),\n  ' ...
                 'end'],State.time{idx_T(i)}, ...
                 State.time{idx_T(i)},State.time{idx_T(i)});
    fprintf(fid,'\n  %s.con(idx,%d) = 0.001;', ...
            c_name,pidx(strmatch(State.time{idx_T(i)}, ...
                                 Param.name,'exact')));
  end
  fprintf(fid,'\n\n  %%variable initialization');
  for i=1:State.n,
    fprintf(fid,'\n  DAE.x(%s.%s) = %s;',c_name,State.name{i},state_init{i});
    fprintf(fid,'\n  %s = DAE.x(%s.%s);',State.name{i},c_name,State.name{i});
  end
  for i=1:nidx,
    fprintf(fid,'\n  %s.dat(:,%d) = %s;',Initl.n+i,c_name,vectorize(Servc.eq{i}));
    fprintf(fid,'\n  %s = %s.dat(:,%d);',Servc.eqidx{iidx(i)},c_name,Initl.n+i);
  end
  for i=1:Initl.n
    fprintf(fid,'\n  %s.dat(:,%d) = %s;',c_name,i, ...
            strrep(Initl.name{i},'_0',''));
  end
  fprintf(fid,'\n\n  %%check limits');
  for i=1:n_xmax
    fprintf(fid,['\n  idx = find(%s > %s_max); if idx, ', ...
                 Comp.name,'warn(idx, '' State variable %s ', ...
                 'is over its maximum limit.''), end'], ...
            State.name{x_max(i)},State.name{x_max(i)}, ...
            State.name{x_max(i)});
  end
  for i=1:n_xmin
    fprintf(fid,['\n  idx = find(%s < %s_min); if idx, ', ...
                 Comp.name,'warn(idx, '' State variable %s ', ...
                 'is under its minimum limit.''), end'], ...
            State.name{x_min(i)},State.name{x_min(i)}, ...
            State.name{x_min(i)});
  end
  for i=1:n_smax
    fprintf(fid,['\n  idx = find(%s > %s_max); if idx, ', ...
                 Comp.name,'warn(idx, '' State variable %s ', ...
                 'is over its maximum limit.''), end'], ...
            Servc.name{s_max(i)},Servc.name{s_max(i)}, ...
            Servc.name{s_max(i)});
  end
  for i=1:n_smin
    fprintf(fid,['\n  idx = find(%s < %s_min); if idx, ', ...
                 Comp.name,'warn(idx, '' State variable %s ', ...
                 'is under its minimum limit.''), end'], ...
            Servc.name{s_min(i)},Servc.name{s_min(i)}, ...
            Servc.name{s_min(i)});
  end
  fprintf(fid,['\n  fm_disp(''Initialization of ',c_name, ...
               'components completed.'')\n']);
end

% ********************************************************************************
% algebraic equations
if ~isempty(algeb_eq)
  if Comp.init, fprintf(fid, '\n case 1 %% algebraic equations\n');
  else,         fprintf(fid, '\n\nswitch flag\n case 1 %% algebraic equations\n');
  end
end

aidx = [1:Algeb.neq];
idx = strmatch('null',Algeb.eq);
aidx(idx) = [];
idx = strmatch('0',Algeb.eq);
aidx(idx) = [];
for i = 1:length(aidx)
  if Buses.n == 1
    a1 = '';
  else
    a1 = num2str(ceil(aidx(i)/2));
  end
  if rem(aidx(i),2)
    fprintf(fid,'\n  DAE.gp = DAE.gp + sparse(%s.bus%s,1,%s,Bus.n,1);', ...
            c_name,a1,vectorize(Algeb.eq{aidx(i)}));
  else
    fprintf(fid,'\n  DAE.gq = DAE.gq + sparse(%s.bus%s,1,%s,Bus.n,1);', ...
            c_name,a1,vectorize(Algeb.eq{aidx(i)}));
  end
end

% ********************************************************************************
% algebraic Jacobians

% substitution of inner service variables
for j = 1:5
  for i = 1:Servc.neq
    if strcmp(Servc.type{i},'Inner') & ~strcmp(Servc.eq{i},'null')
      state_eq = strrep(state_eq,Servc.eqidx{i},['(',Servc.eq{i},')']);
      algeb_eq = strrep(algeb_eq,Servc.eqidx{i},['(',Servc.eq{i},')']);
      servc_eq = strrep(servc_eq,Servc.eqidx{i},['(',Servc.eq{i},')']);
    end
  end
end

if ~isempty(algeb_eq)
  fprintf(fid, '\n\n case 2 %% algebraic Jacobians\n');
end
eqformat = '\n  DAE.J%d%d = DAE.J%d%d + sparse(%s.bus%s,%s.bus%s,%s,Bus.n,Bus.n);';

for j = 1:length(aidx)
  i = aidx(j);
  a1 = 2-rem(i,2);
  if Buses.n == 1
    a2 = '';
  else
    a2 = num2str(ceil(i/2));
  end
  for h = 1:Algeb.n
    type = Algeb.name{h,1};
    if strcmp(type(1), 'V');
      a3 = 2;
      if Buses.n == 1
        a4 = '';
      else
        a4 = type(2:length(type));
      end
    elseif strcmp(type(1:5), 'theta');
      a3 = 1;
      if Buses.n == 1
        a4 = '';
      else
        a4 = type(6:length(type));
      end
    end
    if ~strcmp(char(Gy(i,h)),'0')
      fprintf(fid,eqformat,a1,a3,a1,a3,c_name,a2,c_name,a4, ...
              vectorize(char(Gy(i,h))));
    end
  end
end

% check limits in case of state variable dependancies
Temp = 0;
for i = 1:Servc.neq; Temp = ~strcmp(Servc.limit{i},'None'); break; end
S = 0;
if Temp
  for i = 1:Servc.neq; S = ~strcmp(Servc.type{i},'Input'); break; end
end
if S
  fprintf(fid,'\n');
  for i = 1:length(Servc.eqidx)
    s_var = Servc.eqidx{i};
    for k = 1:Servc.neq; if strcmp(s_var,Servc.name{k}); break; end; end
    if ~strcmp(Servc.type{k},'Input') & ~isempty(findstr(algeb_eq,s_var))
      a = strcmp(Servc.limit{k,1},'None');
      b = strcmp(Servc.limit{k,2},'None');
      if ~a | ~b
        fprintf(fid, ['\n  if (']);
        if ~a
          fprintf(fid,[Servc.name{k},'(i) <= ',Servc.name{k},'_max(i)']);
        else
          fprintf(fid,'(');
        end
        if ~a & ~b
          fprintf(fid,' | '); end
        if ~b
          fprintf(fid,[Servc.name{k},'(i) >= ',Servc.name{k},'_min(i))']);
        else
          fprintf(fid,')');
        end
        fprintf(fid,'\n  end');
      end
    end
  end
end

% ********************************************************************************
% differential & service equations
if ~isempty(state_eq)
  if Comp.init | ~isempty(algeb_eq)
    fprintf(fid, '\n\n case 3 %% differential equations\n');
  else
    fprintf(fid, '\n\nswitch flag\n case 3 %% differential equations\n');
  end
end

for i = 1:Servc.neq
  Temp = Servc.type{i};
  if strcmp(Temp,'Inner')
    s_eq = vectorize(Servc.eq{i});
    fprintf(fid,['\n  ',Servc.name{i},' = ',s_eq,';']);
    if ~strcmp(Servc.limit{i,1},'None')
      fprintf(fid, ['\n  ',Servc.name{i}, ...
                    ' = min(',Servc.name{i},',',Servc.name{i},'_max);']);
    end
    if ~strcmp(Servc.limit{i,2},'None')
      fprintf(fid, ['\n  ',Servc.name{i}, ...
                    ' = max(',Servc.name{i},',',Servc.name{i},'_min);']);
    end
   end
end

for i = 1:State.n
  if strcmp(State.nodyn{i},'Yes')
    fprintf(fid, ['\n  no_dyn_',State.name{i},' = find(',State.time{i},' == 0);']);
    fprintf(fid, ['\n  ', State.time{i}, '(no_dyn_',State.name{i},') = 1;']);
  end
  if strcmp(State.time{i},'None')
    s_eq = vectorize(State.eq{i});
  else
    s_eq = vectorize(['(',State.eq{i},')/',State.time{i}]);
  end
  fprintf(fid, ['\n  DAE.f(',c_name,'.',State.name{i},') = ',s_eq,';']);

  if strcmp(State.nodyn{i},'Yes')
    fprintf(fid, ['\n  DAE.f(',c_name,'.',State.name{i},'(no_dyn_',State.name{i},')) = 0;']);
  end
end

if State.n > 0; if strcmp(State.nodyn{State.n},'Yes'); fprintf(fid, '\n'); end; end

% set hard limits
fprintf(fid,'\n  %% non-windoup limits');
limfor1 = '\n  idx = find(%s >= %s_max & DAE.f(%s) > 0);';
limfor2 = '\n  if idx, DAE.f(%s(idx)) = 0; end';
limfor3 = '\n  DAE.x(%s) = min(%s,%s_max);';
limfor4 = '\n  idx = find(%s <= %s_min & DAE.f(%s) < 0);';
limfor5 = '\n  DAE.x(%s) = max(%s,%s_min);';
for i = 1:State.n
  varidx = [c_name,'.',State.name{i}];
  a = strcmp(State.limit{i,1},'None');
  if ~a
    fprintf(fid,limfor1,State.name{i},State.name{i},varidx);
    fprintf(fid,limfor2,State.name{i});
    fprintf(fid,limfor3,varidx,State.name{i},State.name{i});
  end
  b = strcmp(State.limit{i,2},'None');
  if ~b
    fprintf(fid,limfor4,State.name{i},State.name{i},varidx);
    fprintf(fid,limfor2,State.name{i});
    fprintf(fid,limfor5,varidx,State.name{i},State.name{i});
  end
end

fprintf(fid, '\n');

numdata = Initl.n;
for i = 1:Servc.neq
  Temp = Servc.type{i};
  if okdata & strcmp(Temp,'Inner')
    numdata = numdata + 1;
    fprintf(fid,['\n  ',c_name,'.dat(:,',int2str(numdata),') = ', Servc.name{i},';']);
  elseif strcmp(Temp,'Output')
    numdata = numdata + 1;
    s_eq = vectorize(Servc.eq{i});
    TempT = [c_name,'.dat(:,',int2str(numdata),')'];
    fprintf(fid,['\n  ',TempT,' = ',s_eq,';']);
    zz = ['z(',Servc.name{i},'_',Comp.name,'_idx)'];
    if ~strcmp(Servc.limit{i,1},'None')
      fprintf(fid, ['\n  ',TempT,' = min(',TempT,',',Servc.name{i},'_max);']);
    end
    if ~strcmp(Servc.limit{i,2},'None')
      fprintf(fid, ['\n  ',TempT,' = max(',TempT,',',Servc.name{i},'_min);']);
    end
    fprintf(fid,['\n  ',zz,' = ',zz,' + ',TempT,';']);
  end
end

fprintf(fid, '\n');

% ********************************************************************************
% state variable Jacobians
if ~isempty(state_eq)
  fprintf(fid, '\n\n case 4 %% state variable Jacobians\n');
end

% DAE.Fx
for j = 1:State.n
  if strcmp(State.nodyn{j},'Yes')
    fprintf(fid, ['\n  no_dyn_',State.name{j},' = find(',State.time{j},' == 0);']);
    fprintf(fid, ['\n  ', State.time{j}, '(no_dyn_',State.name{j},') = 1;']);
  end
end
fprintf(fid, '\n');
if State.n, fprintf(fid,'\n  %% DAE.Fx'); end
fxformat = '\n  DAE.Fx = DAE.Fx + sparse(%s,%s,%s,DAE.n,DAE.n);';
for j = 1:State.n
  x_idx1 = [c_name,'.',State.name{j}];
  for i = 1:State.n
    x_idx2 = [c_name,'.',State.name{i}];
    if strcmp(State.time{j},'None')  & ~strcmp(char(Fx(j,i)),'0')
      fxexp = vectorize(char(Fx(j,i)));
    else
      fxexp = ['(',vectorize(char(Fx(j,i))),')./',State.time{j}];
    end
    if ~strcmp(fxexp,['(0)./',State.time{j}])
      fprintf(fid,fxformat,x_idx1,x_idx2,fxexp);
    end
  end
end

fprintf(fid,'\n');

% DAE.Fy
if State.n & Algeb.n, fprintf(fid,'\n  %% DAE.Fy'); end
fyformat = '\n  DAE.Fy = DAE.Fy + sparse(%s,%s,%s,DAE.n,2*Bus.n);';
for j = 1:State.n
  x_idx1 = [c_name,'.',State.name{j}];
  for i = 1:Algeb.n
    type = Algeb.name{i};
    if strcmp(type(1),'V')
      if Buses.n == 1
        x_idx2 = [c_name,'.bus','','+Bus.n'];
      else
        x_idx2 = [c_name,'.bus',type(2:length(type)),'+Bus.n'];
      end
    elseif strcmp(type(1:5),'theta')
      if Buses.n == 1
        x_idx2 = [c_name,'.bus',''];
      else
        x_idx2 = [c_name,'.bus',type(6:length(type))];
      end
    end
    if strcmp(State.time{j},'None') & ~strcmp(char(Fy(j,i)),'0')
      fyexp = vectorize(char(Fy(j,i)));
    else
      fyexp = ['(',vectorize(char(Fy(j,i))),')./',State.time{j}];
    end
    if ~strcmp(fyexp,['(0)./',State.time{j}])
      fprintf(fid,fyformat,x_idx1,x_idx2,fyexp);
    end
  end
end

fprintf(fid,'\n');

% DAE.Gx
if State.n & Algeb.n, fprintf(fid,'\n  %% DAE.Gx'); end
gxformat = '\n  DAE.Gx = DAE.Gx + sparse(%s,%s,%s,2*Bus.n,DAE.n);';
for j = 1:Algeb.neq
  if ~strcmp(Algeb.eq{1},'null')
    type = Algeb.eqidx{j,1};
    if strcmp(type(1),'P')
      if Buses.n == 1
        a_idx = [c_name,'.bus',''];
      else
        a_idx = [c_name,'.bus',type(2:length(type))];
      end
    elseif strcmp(type(1),'Q')
      if Buses.n == 1
        a_idx = [c_name,'.bus','','+Bus.n'];
      else
        a_idx = [c_name,'.bus',type(2:length(type)),'+Bus.n'];
      end
    end
    for h = 1:State.n
      x_idx = [c_name,'.',State.name{h}];
      algexp = vectorize(char(Gx(j,h)));
      if ~strcmp(algexp,'0')
        fprintf(fid,gxformat,a_idx,x_idx,algexp);
      end
    end
  end
end

%if State.n > 0, fprintf(fid, ['\n\n  end']); end

% ********************************************************************************
% non-windup limiters
if n_xmax | n_xmin
  fprintf(fid, '\n\n case 5 %% non-windup limiters\n');
  for i = 1:State.n
    M = ~strcmp(State.limit{i,1},'None');
    m = ~strcmp(State.limit{i,2},'None');
    if M | m
      fprintf(fid, ['\n  idx = find((']);
      if M, fprintf(fid,'%s >= %s_max',State.name{i},State.name{i}); end
      if M & m; fprintf(fid,' | '); end
      if m, fprintf(fid,'%s <= %s_min',State.name{i},State.name{i}); end
      fprintf(fid,[') & DAE.f(',c_name,'.%s) == 0);'],State.name{i});
      fprintf(fid, '\n  if ~isempty(idx)');
      fprintf(fid,['\n    k = ',c_name,'.%s(idx);'],State.name{i});
      fprintf(fid,['\n    DAE.tn(k) = 0;']);
      fprintf(fid,['\n    DAE.Ac(:,k) = 0;']);
      fprintf(fid,['\n    DAE.Ac(k,:) = 0;']);
      fprintf(fid,['\n    DAE.Ac(k,k) = speye(length(idx));']);
      fprintf(fid,['\n  end']);
    end
  end
end

fprintf(fid, '\n\nend\n');

% ********************************************************************************
% warning message function
fprintf(fid,'\n\n%% -------------------------------------------------------------------');
fprintf(fid,'\n%% function for creating warning messages');
fprintf(fid,['\nfunction ',Comp.name,'warn(idx, msg)']);
%fprintf(fid,['\nglobal ',c_name]);
fprintf(fid,['\nfm_disp(strcat(''Warning: ',upper(Comp.name),' #'',int2str(idx),msg))']);

% close component file and return
fclose(fid);
fm_choice(['Function "fm_',Comp.name,'" built.'],2)

% ********************************************************************************
function vect = varvect(vect,sep)

n = length(sep)-1;
if iscell(vect)
  vect = strcat(vect,'#');
  vect = strrep([vect{:}],'#',sep);
  vect(end-n:end) = [];
end