В статье рассматриваются последние достижения в разработках вычислительных процедур моделирования нелинейных радиосхем. В связи с практической актуальностью и высокой вычислительной сложностью рассматриваемая задача является одной из центральных и быстроразвивающихся направлений исследований автоматизации схемотехнического проектирования.
Введение
Последние годы рынок портативных персональных систем радиосвязи резко расширился. Широкая номенклатура проектируемых радиочастотных интегральных схем требует применения средств автоматизации на основных стадиях разработки. Необходимыми являются эффективное и точное моделирование частотных характеристик, искажений, шумовой анализ таких схем, как смесители, фильтры на переключательных конденсаторах, узкополосные усилители.
В отличие от других классов линейных схем, для радиочастотных базовым видом моделирования является нелинейный анализ. Следует сразу отметить, что применение стандартных хорошо развитых средств схемотехнического моделирования, включенных, например, в распространенные программы типа SPICE [1], не позволяет в этом случае, как правило, достигнуть желаемого результата в связи с резко возросшей вычислительной сложностью.
К таким задачам высокой сложности можно отнести, например, необходимую для радиочастотных схем, задачу определения установившихся периодических режимов в нелинейных цепях, задачу моделирования схем с сильно различающимися частотами входных сигналов. В этих задачах требуется моделирование многопериодных нелинейных колебаний на отрезке интегрирования. Практические задачи расчета даже относительно простых цепей, приводят в этом случае к значительным объемам вычислений.
Для оценки возможных типовых затрат моделирования установившихся периодических режимов в работе [2] приведен пример схемы однотранзисторного кварцевого генератора. Схема является высокодобротной с коэффициентом Q=25000. Расчет установившегося режима с помощью программы типа SPICE требует примерно 250000 циклов. Заданная погрешность моделирования в 1% требует выполнения не менее 400 шагов интегрирования в цикле. В результате необходимы расчеты для 10000000 временных точек. По приводимой оценке для ПК PC/486 моделирование периодических режимов для такой схемы занимает 17 часов. При повышенной точности моделирования до 0,1% оценка времени моделирования достигает 30 дней (!)
Таким образом очевиден вывод о необходимости разработки специализированных средств моделирования нелинейных многопериодных цепей, обеспечивающих существенное ускорение в сравнении со стандартными методами моделирования.
Разработка вычислительных схем, позволяющих ускорить моделирование рассматриваемого класса схем имеет уже значительную историю (см., например, [5-12]).
Некоторые результаты нашли применение в промышленно применяемых симуляторах, ориентированных на анализ характеристик радиосхем, например, SPECTRE[1], MDS[4], LIBRA[3] и других. Ниже рассмотрены современные вычислительные схемы для моделирования установившихся периодических режимов в нелинейных цепях (steady-state analysis). Обсуждаются два альтернативных направления:
Преимуществом первого подхода является универсальность, допускающая, в принципе, высокую степень нелинейности элементов цепи. Второй подход получил в последнее время большое распространение, так как очень эффективен для схем, работающих в относительно слабом нелинейном режиме, что характерно для радиочастотных схем, для многовходовых схем с режимом большого сигнала на одном входе. Однако с ростом нелинейности элементов размерность решаемых задач растет очень быстро, что ограничивает возможность применения метода гармонического баланса, даже для практически важных случаев.
Ниже указаны вычислительные схемы, представляющиеся на сегодняшний день наиболее перспективными для решения рассматриваемых задач и позволяющие значительно повысить практические возможности разрабатываемых средств моделирования. Новые подходы, благодаря существенному ускорению, позволили увеличить сложность анализируемых схем.
Методы моделирования во временной области
Будем в дальнейшем использовать два представления математической модели схемы для расчета установившихся периодических режимов. Первый тип модели задается совместной системой
f(x, x, t) = 0, (1a)
x(0) - x(T) = 0. (1b)
Здесь f - система дифференциальных уравнений; x(t) - искомая вектор-функция переменных, определяющих состояние схемы; T - период.
Таким образом определена нелинейная краевая задача, которая может быть решена, вообще говоря, известными методами [13]. Второй тип модели в большей степени отражает свойства электрической цепи [9]:
(2a), (2b)
где u(t) - вектор входных сигналов; v(t) - вектор узловых напряжений; i(v(t)) и q(v(t)) - вектора узловых токов и узловых зарядов (потоков) соответственно.
Этот тип модели позволяет избежать трудностей, связанных с проблемой консервативности зарядов [1].
Применение конечно-разностных методов решения краевой задачи (1), (2) приводит к системе нелинейных уравнений, число неизвестных которой определяется числом точек дискретизации. Например, если обратный метод Эйлера выбран для аппроксимации производной в (1b), мы имеем систему уравнений вида:
(3)
где hj = tj - tj-1, M - число точек дискретизации.
Комбинируя эту систему с условиями периодичности (1b), получаем необходимую для вычисления x(t) и x(0) систему нелинейных уравнений. Применение ньютоновского итерационного процесса приводит к решению на шаге итерации линейной системы вида:
(L + B) (xk+1 - xk ) = - f(xk), (4)
где k - индекс ньютоновской итерации.
Матрицы L и B имеют структуру [20]:
где 
-
емкости и 
- проводимости цепи.
Другим распространненым методом решения задачи (1) или (2) является метод сведения краевой задачи к задаче Коши (shooting method). Этот метод приводит к решению системы нелинейных уравнений (1b) относительно вектора x(0), причем на каждом шаге итерационного процесса вектор x(T) находится решением задачи Коши с соответствующими начальными условиями. Такой подход хорошо согласуется с программами схемотехнического моделирования, поскольку позволяет использовать хорошо отработанные для задач электрического моделирования стандартные методы интегрирования.
Можно отметить, что если систему нелинейных уравнений
x(0) - x(T, x(0)) = 0 (5)
решать методом простой итерации
x(0)k+1 = x(T, x(0)k) (6)
до выполнения условия окончания итерационного процесса по заданной погрешности e:
|| x(0)k+1 - x(0)k|| < e, (7)
то мы, фактически, получаем применение стандартных подходов для поиска периодических решений последовательным интегрированием на большом временном интервале. Основная цель, естественно, сократить исходные затраты моделирования. Для этого, в частности, целесообразно рассмотреть возможность применения метода Ньютона для решения системы (5), требующего меньшего количества шагов. Использование ньютоновского процесса приводит, на шаге итерации, к системе линейных уравнений вида:
где I - единичная матрица.
Основная возникающая проблема в этом случае касается затрат формирования якобиана для (8)
или матрицы чувствительности 
по начальным условиям. Представление этой матрицы в виде
разностей требует многократного решения задачи Коши и практически неэффективно. По этой причине в работах [5, 6] был предложен экономичный
способ вычисления производных . В частности, в работе [6] предлагалось эту матрицу
находить из решения систем уравнений в вариациях:
(9)
Основная идея этого предложения совместить операции, выполняемые при интегрировании системы (1a) и (9), неявными методами. Это приводит к решению на каждом шаге интегрирования системы линейных уравнений с одинаковой матрицей коэффициентов, что позволяет значительно снизить затраты. Алгоритмы такого типа были длительное время базовыми для программ моделирования установившихся колебаний.
Тем не менее, несмотря на экономию операций, методы сведения к задаче Коши с ньютоновской итерацией не удовлетворяют в полной мере возрастающим требованиям к сложностным характеристикам. Их основными ограничениями являются:
Основные направления улучшения алгоритмов связаны с попытками избежать вычисления матрицы , а также сокращения затрат на решение системы с плотными матрицами. Ниже рассматривается развитие этих вычислительных процедур, базирующихся на применении методов подпространств Крылова.
Укажем также предварительно популярный в силу простоты реализации экстраполяционный алгоритм Скилбо [8]. Его идея использовать экстрополяционную схему для ускорения сходимости последовательности, полученной простой итерацией (6):
x(0)0, x(0)1, ... .
В общем случае вычислительную схему можно представить в виде [8]:
x(0) = yn,
x(0)r+1 = X(x(0)r, T), r = 1, 2, ..., m, (10)
yn+1 = EXTRAP (x(0)0, x(0)1, ..., x(0)m).
Преимущества алгоритма заключаются в том, что он позволяет избежать вычисления производных при решении нелинейных задач, позволяет достичь квадратичную сходимость, обеспечивает простоту реализации на базе применяемых стандартных методов интегрирования.
Метод гармонического баланса
В настоящее время метод гармонического баланса является широкораспространенным средством моделирования нелинейных радиосхем [9-11]. Во многих практических случаях он обеспечивает относительно быстрый расчет установившихся режимов как для одночастотных, так и для многочастотных входных воздействий [10].
Этот метод является частотным методом. Его можно определить как метод решения краевой задачи (1) с тригонометрическими базисными функциями. Решение находится в терминах коэффициентов Фурье.
Основные концепции метода гармонического баланса [9]:
Модель, формируемая в частотной области, соответствует исходным временным уравнениям, заданным, например, системой (2). Как правило, исходные уравнения включают также возможность описания распределенных систем и имеют вид [9]:
(11)
После преобразования Фурье эта система может быть переписана в частотной области:
R(V) = I(V) + W Q(V) + Y V + U = 0. (12)
Компоненты этого уравнения являются результатами преобразования F соответствующих временных характеристик:
Fv = V, Fi(V) = I(V), Fq(V) = Q(V), Fu = U, Fy = Y.
Коэффициенты W определены следующим образом [12].
Пусть W = Wmn. Здесь m и n - индексы узлов: m, n = 1, 2, ..., N.
Пусть k и l - индексы частот. Тогда
, где
= W.
Таким образом для формирования модели схемы в виде системы (12) частотные функции I и Q вычисляются для нелинейных элементов схемы в результате следующего цикла:
Основной проблемой метода гармонического баланса является выбор способа численного решения системы нелинейных уравнений(12).
Наиболее распространенные численные алгоритмы базируются на методе Ньютона. В этом случае на итерационном шаге решается система линейных уравнений:
JR (Vj) (Vj+1 - Vj) = -R(Vj) (13)
Здесь JR - гармонический якобиан, который определяется следующим образом:
Обратим внимание на размерности системы (13) и, соответственно, якобиана (14). Она определяется величиной (N*2K), где N - число переменных состояния схемы, K - число гармоник.
При существенной нелинейности анализируемой цепи число K приходится выбирать большим. Соответственно размерность вычислительной модели (12) может оказаться очень высокой. Допустимая размерность решаемой системы (12) является фактическим ограничением метода гармонического баланса для применения к расчету произвольных цепей.
Могут быть предложены самые разнообразные вычислительные схемы ускоренного расчета линейной системы (13), базирующиеся как на прямых методах решения, так и на итеративных. Следует отметить, что применение простейших итеративных методов [16] влечет за собой дополнительную проблему сходимости линейной системы (13), так как эта система, являясь в общем случае моделью активной электрической цепи, не обладает, как правило, желаемыми свойствами диагонального преобладания, положительной определенности или симметричности. Исследования этой проблемы в течение ряда лет (см., например, [9-10]) не привели к результату, который бы более или менее надежно удовлетворял сегодняшним практическим потребностям.
Проблема заключается в радикальном повышении производительности моделирования при естественном желании сохранения надежности сходимости. Лишь результаты недавних исследований позволили сделать шаг в этом направлении.
Проекционные методы на базе подпространств Крылова
Основная конструктивная идея, подтвердившая практическую возможность значительного сокращения вычислительных затрат без ухудшения других численных характеристик, связана с использованием итеративных алгоритмов на базе методов подпространств Крылова. Например, такой подход позволил расширить размер схемы, анализируемый методом гармонического баланса, на 2-3 порядка [17].
Последовательностью Крылова, порожденной матрицей А и вектором Х, называется последовательность векторов x, Ax, A2x, ..., Amx [19]. Соответственно, подпространствами Крылова называют линейные оболочки подсистем последовательности Крылова, т. е. линейные подпространства вида:
Km(x) = span {x, Ax, A2x, ..., Am-1x}.
Проекционные методы решения линейных систем, построенные на базе подпространств Крылова, представляют собой перспективные направления развития современной вычислительной математики. Следует отметить, что основной трудоемкой операцией для них является матрично-векторное умножение (A*x).
Собственно вычислительный алгоритм включает две стадии [19, 23]:
Практически применяемые алгоритмы представляют собой комбинацию методов подпространств Крылова и предобусловливания, которое является вспомогательной операцией, направленной на улучшение вычислительных свойств исходной матрицы без искажения конечного результата. От успешного выбора предобусловливания зависит вычислительная эффективность решения.
С точки зрения применения в схемном моделировании и в анализе многопериодных процессов, в частности, эти методы обладают следующими важными свойствами:
На последнее свойство следует обратить специальное внимание. Исторически подобные алгоритмы применялись, в первую очередь, для симметричных систем, и лишь относительно недавние результаты показали возможность расширения сферы их применения на несимметричные задачи, что является принципиальным вопросом при решении задач схемотехнического моделирования.
Две вычислительные процедуры получили наибольшее распространение в классе обсуждаемых методов [19]: процедура Арнольди и процедура Ланцоша.
На базе процедуры Арнольди разработан для решения линейной системы уравнений A•x = b алгоритм GMRES [22], успешно зарекомендовавший себя в задачах разного типа.
В последние годы целая серия исследований была направлена на применение техники Ланцоша для несимметричных задач. Был получен ряд практических алгоритмов. Среди них можно выделить QMR-алгоритм, разработанный до высокой степени детализации [23].
GMRES и QMR-алгоритмы отличаются разным числом ортонормированных векторов и многими деталями реализации. С точки зрения рассматриваемой задачи и общих вычислительных схем steady-state анализа важно, что для обоих подходов базовой операцией является умножение матрицы на вектор.
Приведем типовую последовательность операций для GMRES-алгоритма при решении системы вида Ах = b [18]:
Задание начального условия х0.
Выбор направления р0 = b - Ax0.
Установить к = 1.
Cледует отметить, что помимо ожидаемой вычислительной экономии при решении высокоразмерных линейных систем, применение обсуждаемых алгоритмов может быть полезно использовано и для модификации собственно вычислительных схем моделирования.
Ниже рассмотрены возможные способы применения для shooting методов и метода гармонического баланса.
В работе [16] предложен подход, повышающий вычислительную эффективность shooting методов при расчете периодических режимов. Его цель - избежать на шаге ньютоновской итерации решения системы уравнений (8) с плотной матрицей стандартными
методами гауссова исключения. Эти методы, как известно, требуют затрат,
пропорциональных ~N3. Матрица 
содержит
N2 элементов. Стоимость одной итерации
при решении линейных систем определяется в основном умножением матрицы
на вектор и оценивается как N2 операций. Если число операций значительно меньше N, то можно расчитывать на результирующий выйгрыш уже в таком варианте.
В работе [18] рассмотрен подход, позволяющий исключить явное формирование якобиана благодаря итеративному процессу, базирующемуся на GMRES-технике. Учитывается тот факт, что в этом случае необходимым для дальнейших вычислений является не сам якобиан 
,
а результат умножения Apk-1. Предлагается рассмотреть разностную форму такого произведения:
(15)
Следовательно, принципиально возможным является определение вектора Apk-1 после интегрирования и вычисления x(Т, x(0) + dpk-1).
В работе [18] предлагается тем не менее использовать уравнение типа (9),
но последовательно, для отдельных направлений, задаваемых вектором pk-1.
При этом вычисляется непосредственно произведение 
в виде вектор-столбцов и операции совмещаются с GMRES-алгоритмом. Как правило, количество итераций существенно меньше N, это и позволяет рассчитывать на достижение экономичности и по приведенным оценкам достичь на практике сложность вычислений порядка О(N) [18]. Выполненные экспериментальные результаты показали, что в среднем ускорение для схемы порядка 400 узлов превышает 10 раз.
Применение итеративных алгоритмов на базе подпространств Крылова для решения линейных систем (13) высокой размерности, формируемых методом гармонического баланса, позволило, как указывалось, радикально повысить возможности этого метода. Действительно, в большинстве практических случаев дополнительной проблемы сходимости для линейных систем не возникает, a число итераций, обеспечивающих заданную погрешность, как правило, много меньше размерности задачи.
В приводимом в работе [17] примере решается система уравнений с 600000(!) неизвестными. Задача решается методом гармонического баланса за 3 часа на вычислительной системе SGI. По приведенной там же оценке, моделирование этого же приемного тракта на тех же вычислительных средствах потребовало бы порядка одной недели процессорного времени.
Следует отметить также, что внедрение нового семейства алгоритмов создает предпосылки, как и во временном анализе, для модификации исходных вычислительных систем с целью дополнительной экономии операций.
Примером может служить изменение порядка операций. Для этого запишем гармонический якобиан (14) в виде [17]:
J = Y + WГСГ-1 + ГGГ-1. (16)
Здесь Г представляет собой оператор преобразования из временной в частотную область, соответственно, Г-1-обратный оператор; матрицы G и C - блочно-диагональные матрицы, вычисленные в точках временной дискретизации (до преобразования Фурье):
При использовании итеративных алгоритмов GMRES или QMR основные затраты уходят на операцию умножения якобиана J на вектор р.
Из выражения
Jp = Yp + WГСГ-1p + ГGГ-1p (17)
видно, что нет необходимости формировать гармонический якобиан J высокой размерности, а затем выполнять умножение всех его элементов на вектор. Более экономичным является предварительное умножение исходных разреженных матриц, вычисленных во временных точках, с последующим применением преобразования Фурье. Таким способом достигается значительная вычислительная экономия. Возможны и другие варианты реализации основных операций.
Быстрое совершенствование численных алгоритмов моделирования нелинейных радиосхем, влияющее на оценку состояния этой области автоматизации электрических расчетов, позволяет сделать лишь условные рекомендации по выбору базовых вычислительных процедур.
1) Несмотря на более высокую универсальность методов моделирования во временной области метод гармонического баланса более практичен для расчета радиочастотных схем, в первую очередь для схем смесителей. Метод гармонического баланса обеспечивает более высокую вычислительную эффективность при моделировании схем с относительно небольшой степенью нелинейности; как правило, к ним относятся радиотехнические схемы.
2) Основные направления увеличения размерностей задач моделирования сложных радиосхем связаны с применением итеративных алгоритмов решения линейных систем на базе подпространств Крылова. Первые практические результаты подтверждают их высокую производительность при сохранении требуемого уровня надежности сходимости.
Заключение
Разработка методов и средств автоматизации схемотехнического моделирования радиотехнических интегральных схем является быстроразвивающимся направлением исследований в области вычислительных методов моделирования. Указанное направление стимулируется практическими потребностями ускоренного развития этой области интегральной схемотехники.
Новые результаты в численном моделировании задач высокой размерности создают предпосылки для радикального увеличения производительности схемных симуляторов при анализе нелинейных радиотехнических схем высокой сложности. Приведенные выше вычислительные процедуры ускоренного расчета нелинейных периодических режимов во временной и частотной областях обеспечивают на практике желаемую точность описания поведения проектируемых схем и являются также вычислительной основой для выполнения таких видов моделирования, как шумовой анализ и анализ искажений.
Литература
1. Kundert K.S. The designer guide to SPICE and SPECTRE. Kluwer Academic Publishers, 1995.
2. Fast and Painless Oscillator Simulation. Avista Design Systems, 1994.
3. Touchstone and Libra for Windows. HP4746 /E4747А Hewlett Packard., 1995.
4. Overview of the RF and Microwave Design System HP851 150B. Hewlett Packard, 1994.
5. Thomas J. Aprille, Timothy N. Trick. Steady-state analysis of nonlinear circuits with periodic inputs. Proceeding of the IEEE, vol. 60, no. 1, pp. 108-114, January 1972.
6. Гурарий М.М., Русаков С.Г., Зарудный Д.И. Моделирование на ЭЦВМ периодических процессов в интегральных схемах. - Автоматика и вычислительная техника, 1973, #1, с. 83-85.
7. Норенков И.П., Евстифеев Ю.А., Маничев В.Б. Ускоренный метод анализа многопериодных схем. - Радиотехника, 1987, #2, с. 71-74.
8. Stig Skelboe. Computation of the periodic steady-state response of nonlinear networks by extrapolation methods. IEEE Transactions on Circuit and Systems, vol. CAS-27, no. 3, pp. 161-175, March 1980.
9. Kundert K.S., White J.K., Sangiovanni-Vincentelli A. Steady-State Methods for Simulating Analog and Microwave Circuit. Kluwer Academic Publishers, 1990.
10. Vittorio Rizolli, Alessandro Lipparini, Alessandra Costanzo. State-of-the-Art Harmonic-Balance Simulation of Forced Nonlinear Microwave Circuit by the Piecewise Technique. IEEE Trans on Microwave Theory and Techniques. Vol. 40, # 1, January 1992.
11. Ланцов В.Н., Меркулов А.С. Алгоритм анализа почти-периодических сигналов в нелинейных радиосистемах. - Радиоэлектроника и системы связи, т.33, 1990, #6, с.12-27.
12. Перминов В.И., Соколов А.Г. Расчет стационарных режимов микросхем. - Микроэлектроника АН СССР, вып. 5, 1984.
13. Холл Д., Уатт Д. Современные методы решения ОДУ. Л., Мир., 1979, с. 312.
14. R. Melville, P.Feldmann, J.Roychowdhury. Efficient multi-tone distortion of analog integrated circuits. Proceedings of the 1995 IEEE Custom Integrated Circuits Conf., May 1995.
15. Бобылев И.А., Бурман Ю.М., Коровин С.К. Оценка погрешности метода гармонического баланса. ДАН, 1991, т 320, # 3, с. 572-576.
16. Ортега Дж. Рейнболдт. Итерационные методы решения нелинейных систем уравнений со многими неизвестными. М., Мир, 1975.
17. P. Feldman, B. Melville, D. Long. Efficient Frequency Domain Analysis of Large Nonlinear Analog Circuits. IEEE 96 Custom Integrated Circuits Conf.
18. R. Telichevesky, K. Kumbert, Jacob K. White. Efficient Steady-State Analysis based on Matrix-Free Krylov-Subspace Methods. Proc. Design Automation Conference, Santa Clara, California, June 1995.
19. Икрамов Х.Д. Разреженные матрицы. В сб. "Итоги науки и техники. Математический анализ", т. 20, 1982, с. 179-260.
20. R. Teleshevesky, K. Kundert, J. White. Fast Simulation Algorithms for RF Circuit. IEEE 96 Custom Integrated Circuits Conf.
21. Saad Y. Krylov subspace methods for solving large unsymmetric linear system. Math. Comp., 37 (1981), pp. 105-126.
22. Saad Y., Shultz. GMRES: A general minimal residual algorithm for solving nonsymmetric linear systems SIAM J. Sci. and Statist. Comput., 1986, 7, n.3, pp. 856-869.
23. Freund R.W., Golub G.H. and Nachtigal N.M. QMR: a quasi-minimal residual method for non-Hermitian linear systems. Numer. Math., 60 (1991), pp. 315-339.