Немного теории

Как было отмечено ранее, при анализе динамических режимов электрических цепей необходимо решить следующую систему уравнений:

,    (1)

где x₁, …, x_n — токи и напряжения в цепи; x'₁, …, x'_m — производные токов и напряжений.

Число производных в системе (1) может быть меньше числа неизвестных токов и напряжений, то есть  m ≤ n, число уравнений в системе всегда равно числу неизвестных x₁, …, x_n. В случае наличия в электрической цепи нелинейных элементов система (1) представляет собой систему нелинейных дифференциально-алгебраических уравнений.

Если определить векторы x=[x₁, x₂, … x_n] ^t и x'=[x'₁, x'₂, … x'_n] ^t, то уравнения (1) можно записать в следующем (матричном) виде:

f(x, x', t) = 0.    (2)

Следует заметить, что если m ≤ n, то в векторе x' будут отсутствовать соответствующие значения.

Если обозначить x(t_j) через x_j, то решение x_n+1 уравнения (2) в точке t = t_n+1 должно удовлетворять уравнению (здесь n — номер точки решения)

f(x_n+1, x'_n+1, t_n+1) = 0.    (3)

В соответствии с формулой дифференцирования назад (BDF) значение производной x'_n+1 аппроксимируется с использованием x_n+1 и K последних значений x_n, x_n-1, …, x_n-K+1

,   (4)

где α_j — коэффициенты BDF; h — текущий шаг интегрирования; K — порядок BDF. В правой части выражения (4) неизвестной является только текущее значение x_n+1, поэтому обозначим как функцию g(x_n+1).

Подставив выражение (4) в (3), получим уравнение

f(x_n+1, g(x_n+1), t_n+1) = 0,

которое является системой нелинейных алгебраических уравнений относительно неизвестной x_n+1. Для решения этой системы может быть применен, например, метод Ньютона–Рафсона. После определения x_n+1 можно аналогичным образом найти x_n+2, x_n+3 и т.д.

В начальный момент времени t = t₀ известны лишь начальные значения x₀ = x(t₀), поэтому в выражении (4) следует взять порядок K = 1 и найти x₁. Затем, зная x₀ и x₁ можно вычислить x₂, применяя порядок K ≤ 2. Поступая таким образом, можно последовательно выйти на требуемый порядок BDF, в большинстве случаев используется K ≤ 6.

Таким образом при наличии в электрической цепи нелинейных элементов расчет динамического режима с использованием BDF можно свести к решению системы нелинейных алгебраических уравнений методом Ньютона-Рафсона на каждом шаге интегрирования. При этом нужно помнить, что итерационный процесс сходится, когда начальная приближенная оценка  достаточно близка к решению x_n+1 уравнения (3). В нашем случае, когда известны последние K + 1 значений предыдущих решений, можно вычислить хорошее начальное приближение с помощью формулы прогноза BDF:

,    (5)

Коэффициенты γ₁, γ₂, …, γ_K+1, входящие в выражение (5), при неравномерном шаге интегрирования должны вычисляться для каждого нового момента времени, формулы расчета коэффициентов приведены в [1] и [2]. При постоянном шаге интегрирования коэффициенты прогноза BDF представляют собой константы и не зависят ни от времени, ни от шага интегрирования, их значения можно найти в [3].

В начальный момент времени t = t₀, когда известно лишь начальное значение x₀, прогнозное значение следует принять .

Расчет параметров дискретной модели нелинейной индуктивности

Рассмотрим особенности расчета параметров дискретной модели нелинейной индуктивности,  характеристика которой задана в виде зависимости Ψ = f(i), где Ψ — потокосцепление катушки. Характеристика может быть задана либо  в табличной форме, либо в виде аналитического выражения (по ссылкам приведены примеры, как это выполнено в LazDiscret2). В случае применения таблично-заданной характеристики для аппроксимации можно использовать, например, кубический сплайн, который позволяет и выполнять интерполяцию характеристики между опорными точками.

Учитывая связь напряжением на индуктивности и потокосцеплением, а также выражение BDF, для момента времени t_n+1 можно записать

.   (6)

Значение потокосцепления на k+1-ом шаге итерации решения методом Ньютона-Рафсона можно найти [4]

,

здесь верхний индекс у переменных обозначает номер итерации решения методом Ньютона-Рафсона, а нижний индекс — момент времени для интегрирования с помощью BDF.

Значения потокосцепления и его производной по току могут быть найдены по заданной характеристике Ψ = f(i):

  и  ,

где calcValue и calcDeriv — функции вычисления интерполяционных значений величины потокосцепления и его производной для заданной величины тока. Для таблично-заданной характеристики эти значения вычисляются с помощью принятой аппроксимирующей функции, в случае аналитически заданной характеристики — по функции характеристики, определенной на этапе редактирования параметров.

Для первого шага итерации в качестве предыдущего значения тока (начальное приближение) берется значение, вычисленное по формуле прогноза

,

для которой необходимо сохранение значений тока на предыдущих K+1 шагах интегрирования (K — порядок BDF).

В итоге для момента времени t_n+1 на k+1 шаге итерации выражение (6) можно записать в следующем виде

,

где

— эквивалентное сопротивление дискретной модели нелинейной индуктивности на k+1 шаге итерации для момента времени t_n+1;

— эквивалентная ЭДС дискретной модели нелинейной индуктивности на k+1 шаге итерации для момента времени t_n+1.

Итеративный процесс для момента времени t_n+1 необходимо повторять до достижения условия сходимости решения. Для исключения зацикливания при расходимости итеративного процесса нужно ограничить максимальное числ шагов итерации.

Аналогичные модели можно получить для нелинейной емкости с заданной характеристикой Q = f(u) и для нелинейной взаимной индуктивности с характеристикой материала сердечника B = f(H) или Φ = f(i_нам), где i_нам — ток намагничивания.

Алгоритм расчета динамического режима нелинейной цепи

Алгоритм расчета нелинейной цепи при постоянном шаге интегрирования упрощенно выглядит следующим образом:

//расчет/установка начальных значений
calculateInitial;    
//первый шаг интегрирования
global_Time := delta;
//начальный порядок
BDF_Order := 1;
BDF_Order_Changed := true;

//цикл расчета во временной области
repeat
  //проверка изменения порядка BDF
  if BDF_Order_Changed then begin
    //расчет эквивалентных параметров линейных элементов
    calculateLinearDynBranchesPars;
    BDF_Order_Changed := false;
  end;

  //расчет текущих alpha сумм
  calculateAlphaSumm;
  //расчет значений прогноза
  calculatePredictor;
  //счетчик итераций = 0
  iterCount := 0;
  //отметка о необходимости установки прогнозного значения
  include(modeSet, setPredictor);

  //цикл расчета методом Ньютона-Рафсона
  repeat
    //формирование матриц ММУП
    matrixSetup;
    //решение системы ММУП
    matrixSolve;
    //расчет токов и напряжений элементов
    calculateCurrent;
    //проверка сходимости
    if convergence_check then break;
    iterCount += 1;  
    exclude(modeSet, setPredictor);
  until iterCount > maxIterCount; //конец цикла итераций

  if iterCount > maxIterCount then error; //превышение макс. числа итераций
  //сохранение переменных состояния
  saveStateVars;
//сохранение управляющих переменных нелинейных элементов
  saveNLControlVars;
  //следующий шаг интегрирования
  global_Time += delta;

  //проверка необходимости увеличения порядка BDF
  if BDF_Order < maxBDF_Order then begin
    BDF_Order += 1;
    BDF_Order_Changed := true;
  end;

until global_Time > end_Time; //конец цикла расчета BDF

Возможная реализация соответствующих процедур для модели нелинейной индуктивности описана далее.

Предположим, что для работы с BDF имеются следующие процедуры:

function BDF.AlphaSumma(const aSavedVar: TStateArray): TFloat;
var i: integer;
begin
  result := 0.0;
  for i := 1 to BDFOrder do
    result += aSavedVar[i] * Alpha[i+1];
end;

function BDF.GammaSumma(const aSavedVar: TStateArray): TFloat;
var i: integer;
begin
  result := 0.0;
  for i := 0 to BDFOrder do
    result += aSavedVar[i] * Gamma[i];
end;

procedure BDF.addToStateArray(var aSavedVar: TStateArray; aValue: TFloat);
var i: integer;
begin
  for i := High(aSavedVar) downto 1 do
    aSavedVar[i] := aSavedVar[i-1];
  aSavedVar[0] := aValue;
end;

Модель нелинейной индуктивности имеет следующие переменные

FStateVar — сохраненные значения Ψ в моменты времени t_n, t_n-1, t_n-2, ...;
FCtrlVar — сохраненные значения i_L в моменты времени t_n, t_n-1, t_n-2, ...;
ibk — ток i_L k-того шага итерации в момент времени t_n;
ipred — прогнозное значение тока i_L;
eeqk — эквивалентная ЭДС модели для k+1-го шага итерации в момент времени t_n;

Также процедуры, используемые на различных этапах расчета:

//calculateAlphaSumm; ------------------------------------------------------------------------------
   eeqk := - BDF.AlphaSumma(FStateVar) / integrStep;

//calculatePredictor; ------------------------------------------------------------------------------
   ipred := BDF.GammaSumma(FCtrlVar);

//matrixSetup; -------------------------------------------------------------------------------------
   if setPredictor in modeSet then
     //первый шаг итерации - используем прогноз
     ibk := ipred;   
   psik  := evalTree.calcValue(ibk); //расчет потокосцепления
   ldifk := evalTree.calcDeriv(ibk); //расчет производной потокосцепления
   reqk  := -(BDF.Alpha0 * ldif) / integrStep; //эквивалентное сопротивление
   eeqk   += -(BDF.Alpha0 * (psi - ldif * ibk)) / integrStep; //эквивалентная ЭДС

   //формирование матрицы проводимостей ММУП
   ymat[pnode, eqnN] += 1.0;    ymat[nnode, eqnN] += -1.0;
   ymat[eqnN, pnode] += 1.0;    ymat[eqnN, nnode] += -1.0;
   ymat[eqnN, eqnN]  -= reqk;
   //формирование ветора правой части ММУП
   rhs[eqnN] += eeqk;

//calculateCurrent; --------------------------------------------------------------------------------
   ubk := solution[pnode] - solution[nnode]; //расчет напряжения на ветви
   ibk := solution[eqnN];  //определение нового тока

//saveStateVars; -----------------------------------------------------------------------------------
    BDF.addToStateArray(FStateVar, evalTree.calcValue(ibk));

//saveNLControlVars; -------------------------------------------------------------------------------
    BDF.addToStateArray(FCtrlVar, ibk);

Аналогично могут быть построены процедуры для моделей других нелинейных элементов.

К началу

Литература

1. Чуа Л.О., Лин Пен-Мин. Машинный анализ электронных схем: Алгоритмы и вычислительные методы. - М.: Энергия, 1980

2. Влах И., Сингхал К. Машинные методы анализа и проектирования электронных схем. - М.: Радио и связь, 1988

3. Нерретер В. Расчет электрических цепей на персональной ЭВМ М.: Энергоатомиздат 1991

4. L. W. Nagel. SPICE2: A computer program to simulate semiconductor circuits. Univ. of California, Berkley, ERL Memo ERL M520(May). 1975