ББК 32.973 И46 УДК 681.14 Ильин В. Н., Коган В. Л. И46 Разработка и применение программ автомати- зации схемотехнического проектирования. — М.: Радио и связь. — 1984. — 368 с., ил. В пер.: 1 р. 30 к. Рассмотрены постановка задач автоматизированного проектирова- ния радиоэлектронных устройств, методика формального описания этих задач, принципы построения моделей элементов и алгоритмы расчета устройств различных типов. Приведены способы разработки структуры программ, их отладки и технология использования. Для инженерно-технических работников, занимающихся разработ- кой систем автоматизированного проектирования РЭА и их исполь- зованием. „2401000000-195 " 046(01 )-84 63-84 ББК 32.973 6Ф7.3 Рецензенты: канд. техн. наук А. Я. Архангельский и канд. техн. наук С. Г. Русаков Редакция литературы по кибернетике и вычислительной технике ВАЛЕРИЙ НИКОЛАЕВИЧ ИЛЬИН ВЛАДИМИР ЛЕОНИДОВИЧ КОГАН РАЗРАБОТКА И ПРИМЕНЕНИЕ ПРОГРАММ АВТОМАТИЗАЦИИ СХЕМОТЕХНИЧЕСКОГО ПРОЕКТИРОВАНИЯ Редактор Н. Г. Давыдова Обложка художника В. В. Волкова Художественный редактор Н. С. Шеин Технический редактор Т. Н. Зыкина Корректор Т. Г. Захарова ИБ № 406 Сдано в на"юр 24.04.84 Подписано в печать" 11.07.84 Т-ОБ458 Формат 84 X 108'/32 Бумага кн.-журнальная Гарнитура литературная Печать высокая Усл. печ. л. 19,3' Усл. кр.-отт. 19,32 Уч.-изд. л. 20,4 Тираж 6000 экз. Изд. № 20143 Зак. № 3469 Цена 1 р. 30 к, Издательство <Радио и связь». 101000 Москва, Почтамт, а/я 693 Ордена Октябрьской Революции и ордена Трудового Красного Знамени Пер- вая Образцовая типография имени А. А. Жданова Союзполиграфпрома при Государственном комитете СССР по делам издательств, полиграфии и книж- ной торговли. 113054, Москва, М-54, Валовая, 28 © Издательство «Радио и связь», 1984 ПРЕДИСЛОВИЕ Процесс проектирования сложных радиоэлектронных изделий в общем случае включает этапы структурного, функционально-логического, схемотехнического, компо- нентного и конструкторского проектирования. Этап схе- мотехнического проектирования играет важную роль в этом процессе и состоит в тщательной проработке прин- ципиальной электрической схемы устройства в целом и отдельных его узлов. Эволюция программного обеспечения схемотехниче- ского проектирования не только существенно расширила возможности автоматизированного расчета схем на ЭВМ, но и одновременно поставила перед разработчи- ками программ и их пользователями ряд сложных про- блем. Для разработчиков программ автоматизации схемо- технического проектирования (АСхП) основную про- блему можно сформулировать следующим образом: как создать современную программу АСхП, максимально удовлетворяющую запросы пользователя и наиболее эффективно использующую программные и аппаратные возможности существующих вычислительных средств. Для решения этой проблемы необходимо разработать: лаконичный и информативный входной язык, удоб- ный для описания как схемы, так и различных заданий на ее проектирование, а также для управления работой программы; гибкую макроструктуру программы (разбиение на отдельные модули), допускающую ее развитие и быструю модификацию, обеспечивающую возможность решения различных задач проектирования и позволяющую ком- поновать с помощью операционной системы высокоэф- фективные в отношении вычислительных затрат овер- лейные структуры; компактные и гибкие микроструктуры отдельный программных модулей, в том числе различных библио- тек, каталогов и др., отличающиеся высокой скоростью реализации заложенных в них алгоритмов; при этом нужно использовать наиболее эффективные приемы про- граммирования. Однако практика показывает, что недостаточно толь- ко разработать программу, необходимо ее внедрить, обеспечить эксплуатацию в каждой конкретной органи- зации или коллективе пользователей. Отсюда вытекает вторая проблема, относящаяся к той группе пользова- телей, которая занимается внедрением и обслуживани- ем уже созданных программ, т. е. пользователей-сопро- водителей: как обеспечить квалифицированную эксплу- атацию программы АСхП, позволяющую удовлетворить запросы различных проектировщиков схем. Для этого сопроводитель программы должен уметь: создавать нужные конкретному проектировщику ре- жимы работы программы и компоновать соответствую- щие им модульные структуры; диагностировать возможные неудачи и сбои при ра- боте программы; разрабатывать отдельные модели элементов, сигна- лов и др. в соответствии с используемой элементной базой и в той форме, которая принята в эксплуатируе- мой программе; оперативно пополнять библиотеки моделей элемен- тов, сигналов и макромоделей функциональных узлов разработанными или известными моделями и макромо- делями. Одно из важных мероприятий по сопровождению программ состоит в разработке и освоении пользовате- лями методики определения параметров имеющихся в программе моделей и макромоделей элементов и узлов схем. Практика показывает, что эффективность использо- вания программ АСхП в значительной степени опреде- ляется квалификацией групп сопровождения, участни- ки которых должны владеть перечисленными навыками. Наконец, окончательный эффект применения про- грамм АСхП зависит от того, насколько правильно они будут использоваться в задачах проектирования схем наиболее многочисленной группой пользователей — не- посредственными проектировщиками схем. Основную проблему, стоящую перед ними, можно сформулировать следующим образом: как спроектировать схему с по- мощью данной программы. Для этого проектировщик схемы должен: уметь описать схему и задание на ее проектирова- ние на входном языке программы; знать функциональные возможности программы, т.е. какие конкретно задачи проектирования в принципе можно решать на ЭВМ с ее помощью; знать информационное обеспечение программы, т.е. состав библиотек моделей элементов, сигналов и др. и их параметры; иметь представление о затратах машинного времени и памяти, необходимых для решения различных задач проектирования, во избежание постановки нереальных задач слишком большого размера или требующих для своего решения непомерно большого времени; владеть основными приемами настройки программы на задачу и задачи на программу для получения тре- буемого результата с минимальными затратами времени и памяти и достаточной точностью. Книга адресована всем трем указанным выше груп- пам лиц—разработчикам программ, сопроводителям и пользователям. Для разработчиков программ наиболь- ший интерес представляют гл. 2, 5—8, в которых при- ведены основные сведения об алгоритмах АСхП и тех- нологии разработки программ, начиная от эскизного проекта и кончая процессом отладки. Главы 1, 3, 4, 9 ориентированы в основном на пользователей и содер- жат необходимые им сведения о терминологии, исполь- 5 зуемой в задачах АСхП, типах программ и их возмож- ностях, моделях типовых элементов схем, входных язы- ках программ АСхП и методологии применения про- грамм. Значительная часть книги представляет интерес и для сопроводителей программ АСхП—это методика составления моделей и макромоделей элементов и изме- рения их параметров (гл. 3), тестирование программ (гл. 8), методология их запуска и сопровождения (гл. 9). Разумеется, трудно провести четкую границу интере- сов каждой из трех указанных групп лиц, занимающих- ся разработкой и использованием программ АСхП. По- этому, например, гл. 2, где излагаются основные алго- ритмы АСхП, написана так, чтобы ее содержание было полезным и интересным не только для разработчика программ, но и пользователя, который мог бы легко составить представление о существующих алгоритмах АСхП, их особенностях и основных проблемах решения типовых задач проектирования на ЭВМ. Следует отметить, что проблема автоматизации про- ектирования стоит сейчас во многих отраслях народно- го хозяйства. Вследствие того, что разнообразные по физической природе проектируемые объекты имеют ча- сто одинаковое математическое описание, например в виде системы обыкновенных дифференциальных уравне- ний, излагаемые в книге алгоритмы и методы разработ- ки программ АСхП можно использовать не только в радиоэлектронике, но и в других областях. Это должно способствовать сокращению сроков разработки про- грамм для САПР в этих областях. В настоящее время в ряде организаций эксплуатиру- ются различные программы АСхП [I], в том числе отраслевого значения («Рапира», САМРИС-2 и др.). Однако многообразие задач АСхП, различие используе- мых для проектирования вычислительных средств, от- сутствие общесоюзных стандартов на модели элемен- тов, языки проектирования и другие причины не позво- ляют ориентировать всех пользователей на одну или даже несколько разработанных программ. В связи с этим книга посвящена не детальному описанию каких- либо конкретных программ, а обобщению накопленного сейчас опыта их разработки и использования. В отличие от ряда известных монографий, посвящен- ных описанию теоретических алгоритмов АСхП, ,в дан- ной книге основное внимание уделено вопросам разра- ботки самих программ, входному языку, моделям ти- повых элементов схем, методике отладки программ и их использования, пока еще недостаточно освещенным в литературе. Учитывая сложившуюся к настоящему времени практически повсеместно тенденцию к приме- нению в программах АСхП в качестве базисного—ме- тода узловых потенциалов и его модификаций, в качест- ве основных численных алгоритмов — методов неявного интегрирования дифференциальных уравнений в сочета- нии с методом Ньютона и учетом разреженности матриц коэффициентов схемы, а в качестве основного языка программирования — языка Фортран, авторы сочли це- лесообразным ориентировать изложение в части разра-' ботки программ именно на эти их особенности. Кроме того, описание всех моделей компонентов схем — бипо- лярных и полевых транзисторов, трансформаторов, рас- пределенных структур, зависимых источников — доведе- но до такого уровня детализации, который позволяет включать эти модели в программы АСхП, основанные на методе узловых потенциалов, без каких-либо серь- езных доработок. В заключение авторы выражают признательность канд. техн. наук Г. П. Мозговому, совместно с кото- рым написан § 9.8, и рецензентам—кандидатам техн. наук А. Я. Архангельскому и С. Г. Русакову, замечания которых способствовали существенному улучшению со- держания книги. ГЛАВА 1 ВВЕДЕНИЕ В этой главе приведены основные необходимые поль- зователю сведения о программах АСхП, их структуре, типах задач АСхП и используемой терминологии, от- ражена степень автоматизации решения различных схе- мотехнических задач для различных классов схем. 1,1. Общая характеристика программ АСхП В отличие от разработчика программ АСхП, для ко- торого важны структура и состав модулей, информаци- онные и управляющие связи программы, пользователя в первую очередь интересуют ее принцип действия и функциональные возможности. С этой точки зрения ра- боту любой программы можно представить следующим образом (рис. 1.1). 1. Описание схемы, включая описание задания на ее проектирование, составленное на входном языке про- граммы (языке проектирования, обычно отличающемся от языка программирования—Фортрана, ПЛ/1 и др.), с помощью подпрограмм ввода вводится в память ЭВМ и запоминается в соответствующих массивах. 2. С помощью подпрограмм, в совокупности назы- ваемых трансляторами входного языка, описание схемы «расшифровывается» и перерабатывается. Целью рас- шифровки является проверка правильности составления описания и выдача на печать сообщений об ошибках (диагностика ошибок), а целью переработки—переком- поновка данных о схеме из формы, удобной для состав- ления ее описания, в новую форму, удобную для работы программы. Эта новая форма называется обычно внут- ренним форматом данных. 3. На основе внутреннего формата составляется ма- тематическая модель схемы (ММС) не в виде исходных уравнений типа Р(л;)==0 или dx/dt^f{x, t) и обычно не в виде формул решения этих уравнений, поскольку, как правило, их невозможно получить, а в виде векто- ров и матриц, входящих в состав реализованного в про- грамме численного алгоритма решения исходных урав- нений (например, матрицы Якоби F' (х) и вектора не- вязок F{x) для алгоритма Ньютона решения исходного уравнения F(x)==0). Векторы и матрицы представляют- Представление информации Обработка информации Исходное описание схемы на входном языке Подпрограммы ввода Исходные массивы данных в памяти ЭВМ Подпрограммы трансляции описания схемы Описание схемы во внутреннем формате Типы и параметры моделей -—————^ Топология схемы Подпрограммы формирования математической модели схемы Библиотека моделэй элементов схем Векторы и матрицы, входящие • алгоритм расчета модели схемы Подпрограммы реализации алгоритма расчета модели схемы Массивы результатов расчета модели схемы Подпрограммы обработки результатов расчета схемы Массивы выходных параметров схемы Подпрограммы анализа Подпрограммы оптимизации Рис. 1.1. АСхП Представление и обработка информации в программе ся в памяти ЭВМ соответствующими массивами число- вых данных (а не формул!) и обычно вычисляются за- ново на каждой новой итерации алгоритма. 4. Эти векторы и матрицы вычисляются программой автоматически путем обращения к библиотекам мате- матических моделей элементов схем, моделей сигналов и др. Процесс составления ММС состоит в последова- тельном обращении к подпрограммам моделирования элементов схемы и определения вклада каждого эле- мента в общую модель схемы. Этот вклад может быть, например, током элемента в составе узлового тока или напряжением на элементе в составе контурного напря- жения и т. д. Описание элементов схемы, входящее в описание схе- мы, участвует в составлении ММС различным образом: с одной стороны, типы и параметры элементов схем используются при обращении к моделям элементов (че- рез каталог библиотеки моделей) и при расчете моде- лей, а с другой — топология схемы (узлы включения элементов в схему) используется при формировании ма- триц и векторов ММС из вкладов отдельных элементов схемы. 5. На основе заложенного в программу алгоритма рассчитывается модель схемы с учетом указанных поль- зователем в задании на расчет директив (время расче- та, точность и т. д.). Результаты расчета (токи, напря- жения, потенциалы) запоминаются в соответствующих массивах и затем обрабатываются с целью вычисления выходных параметров (фронтов, длительностей импуль- сов, параметров частотных характеристик и др.). Рассмотренная программа называется моделирую- щей и образует основную часть более сложных про- грамм анализа и оптимизации, в которых обращение к моделирующей программе как подпрограмме выпол- няется на каждом шаге анализа или оптимизации. Отметим, что описанная последовательность работы программы не совпадает с последовательностью вызо- вов подпрограмм внутри программы. После формирова- ния внутреннего формата вызывается старшая (нижняя на рис. 1.1) из участвующих в работе подпрограмм, на- пример подпрограмма расчета выходных параметров, внутри ее вызывается подчиненная подпрограмма фор- мирования массивов токов или напряжений, необходи- мых для расчета выходных параметров, и т. д. вплоть до самых младших — подпрограмм расчета моделей эле- 10 ментов. Более подробно вложенная структура программ рассматривается в гл. 6. Программные средства АСхП прошли большой путь развития, отражающий совершенствование теории и тех- нических средств. Можно условно выделить три поко- ления программ АСхП, в которых математический ап- парат основан на решении обыкновенных дифференци- альных уравнений (ОДУ). К первому поколению отно- сятся программы, характеризующиеся ограничениями на шаг интегрирования из-за использования явных мето- дов численного решения ОДУ, описывающих переход- ные процессы в схемах. Большинство программ перво- го поколения основано на методе переменных состоя- ния. Ко второму поколению принадлежат программы, не имеющие ограничений на шаг интегрирования, так как в них использовались неявные методы интегрирования и др. В большинстве программ второго поколения учи- тывалась разреженность матриц уравнений схемы, име- лась значительно большая по сравнению с первым по- колением программ свобода в использовании различ- ных математических моделей активных элементов, представляемых, как правило, не схемами замещения из двухполюсных элементов, а многополюсниками. Оба первых поколения программ характеризуются ограни- ченными функциональными возможностями, сводящи- мися в большинстве случаев к расчету статического ре- жима и переходных процессов. Наконец, к третьему поколению можно отнести на- ходящиеся сейчас в эксплуатации сложные комплексы I программ [I], вобравшие в себя все достоинства про- i 01 грамм второго поколения, но отличающиеся более ши- рокими функциональными возможностями, в том числе возможностями решения задач анализа — статистиче- ского, температурного и т. д. и, что особенно важно, за- дач оптимизации параметров схем, решаемых в едином вычислительном процессе с задачами расчета и ана- лиза. Возможны и другие способы деления программ на поколения [2, З], однако в любом случае это деление весьма условно. Одна и та же программа по разным признакам (функциональным возможностям, алгорит- мам, структуре) может быть отнесена к разным поко- лениям. Важно отметить, что совершенствование и усложнение программного обеспечения АСхП во мпо- 11 гом определяется не только развитием теории расчета схем, но и разработкой новой вычислительной техники 1ак, весьма ограниченный объем памяти ЭВМ пео- вого и второго поколения оказал серьезнейшее влияние на алгоритмы программ первого поколения. Основным требованием к этим алгоритмам был экономный расход памяти, что привело, в частности, к широкому исполь- зованию явных методов численного решения ОДУ Появление ЭВМ третьего поколения с большим ооъемом оперативной памяти открыло дорогу примене- нию в программах АСхП новых, более быстрых,^тре- бующих больших затрат памяти неявных методов реше- ния ИДУ, позволило перейти от решения отдельных частных задач расчета схем (статика, динамика) к ре- шению комплексных задач (многовариантный анализ оптимизация). Использование на ЭВМ третьего поколе-' ния современных операционных систем ОС ЕС, ДОС ЕС, ДОС АСВТ привело к появлению мощных комплек- сов программ АСхП для проектирования схем в пакет- ном режиме. С помощью операционных систем можно компоновать различные конфигурации этих комплексов соответствующие поставленным задачам проектировав ния. • Разработка периферийного оборудования ЭВМ (дис- плеев, графопостроителей и др.) вызвала дальнейшее развитие программ АСхП, в которых широко использу- ется дисплейный ввод и графический вывод информа- ции, а также диалоговый режим работы. По современным представлениям каждая САПР в том числе система АСхП, имеет методическое, програм- мное, информационное, техническое и организационное обеспечение (рис. 1.2). Как следует из рис. 1.2, разра- ботка и эксплуатация системы АСхП представляет сложную научную, техническую и организационную за- дачи. 1.2. Возможности автоматизации схемотехнического проектирования различных классов электронных схем В каждой программе АСхП заложен ограниченный круг алгоритмов, рассчитанных на проектирование лишь определенного типа схем. В связи с этим проек- тировщику полезно знать возможности автоматизации схемотехнического проектирования различных электрон- ных схем, классификация которых приведена в табл. 1.1. 12 Таблица 1.1 Тип электронных схем Назначение Вид сигнала Радиочастотные Генераторы, усилители, фильтры, модуляторы, детекторы, смесители, СВЧ-устройства и т. д. Гармонический немодулированный несущий сигнал. Гармонический модулированный несущий сигнал Видеочастотные (с негармониче- Генераторы, усилители, логические схемы, фор- Низкочастотный непрерывный сигнал. Одиночные импульсы. Модулированные" импульсные последовательности Для радиочастотных схем наиболее полно разрабо- таны программы расчета радиочастотных линейных усилителей с сосредоточенными параметрами. Расчет выполняется в частотной области. Результатом расче- та являются амплитудно-частотные и фазо-частотные характеристики. Кроме того, имеются программы ана- лиза устойчивости радиочастотных усилителей. В значительно меньшей степени разработаны про- граммы расчета модуляторов и детекторов, в которых высокочастотный несущий сигнал модулирован по ам- плитуде, частоте или фазе в соответствии с законом изменения низкочастотного информационного сигнала. Анализ этих схем во временной области требует выде- ления соответствующей амплитудной, частотной или фа- зовой огибающей несущего сигнала. Современные про- граммы АСхП позволяют выделять только амплитудную огибающую, поэтому схемы с частотно- и фазомодули- ровапными сигналами рассчитываются не во времен- ной, а в частотной области путем определения по задан- ному спектру входного сигнала спектра выходного сиг- нала. Этот же подход используется для расчета в ча- стотной области схем с амплитудно-модулированным сигналом, а также нелинейных радиочастотных схем. Основные трудности расчета нелинейных радиоча- стотных схем состоят в том, что расчет во временной области требует больших затрат времени, а расчет в частотной области может привести к потере точности, 14 если для экономии времени учитывать в спектрах вход- ного и выходного сигналов небольшое число частот. Наконец, следует отметить, что разработанные про- граммы АСхП радиочастотных схем предназначены в основном для проектирования схем с сосредоточенными параметрами. Имеющиеся программы расчета радиоча- стотных схем с распределенными параметрами носят спе- циализированный характер и требуют дальнейшей уни- версализации и сокращения времени счета. В целом состояние разработок программ для АСхП радиочастотных схем не позволяет пока полностью удов- летворить проектировщиков этого класса схем. Это вы- звано более поздним началом работ в области АСхП ра- диочастотных схем по сравнению, например, с работа- ми по АСхП переключательных схем. Однако достигну- тые результаты и развернутый фронт исследований по- зволяют надеяться на успешное решение задач АСхП радиочастотных схем. Второй обширный класс электронных схем состав- ляют видеочастотные схемы, которые часто называют аналоговыми, но это название нельзя признать удач- ным, поскольку многие «аналоговые» схемы работают в составе чисто цифровых устройств. К тому же трудно согласиться с тем, что триггер является аналоговой схемой. Наиболее полно разработаны программы расчета видеочастотных схем, работающих с низкочастотными непрерывными и с одиночными импульсными сигнала- ми. Обычно линейные усилительные схемы рассчитыва- ются как во временной, так и в частотной области, а остальные типы видеосхем — только во временной. Следует отметить, что качество функционирования видеосхем, работающих с модулированными импульс- ными последовательностями, может быть определено в результате анализа не одного, а некоторой достаточно Длинной последовательности импульсов, играющей роль высокочастотного заполнения. По существу эти схемы являются импульсными системами. Анализ таких си- стем с помощью обычных программ расчета импульс- ных схем во временной области занимает много време- ни, а специальные методы и программы АСхП этих систем разработаны недостаточно. По типу параметров схемы делятся на схемы с рас- пределенными и сосредоточенными параметрами. Схе- мы с распределенными параметрами описываются урав- 15 нениями в частных производных, численное решение ко- торых на ЭВМ требует больших затрат машинного вре- мени. Поэтому всегда, если это возможно, стараются представить схемы с распределенными параметрами ап- проксимирующими схемами с сосредоточенными пара- метрами, расчет которых выполняется значительно бы- стрее. Все существующие схемы можно разделить на три основных вида: 1) на дискретных компонентах—тран- зисторах в корпусном или бескорпусном исполнении; 2) на дискретных функциональных элементах — инте- гральных микросхемах; 3) полупроводниковые инте- гральные схемы, изготавливаемые целиком в едином технологическом процессе. Составление и расчет урав- нений этих трех видов схем может быть выполнен с помощью одних и тех же программ АСхП, однако ин- формационное обеспечение этих программ и маршруты автоматизированного проектирования схем различны (см. гл. 9). По характеру уравнений, описывающих процессы и состояния схем, схемы разделяются на линейные и не- линейные. Наиболее сложен расчет нелинейных радио- частотных схем, методы которого разработаны недоста- точно. Расчет линейных схем выполняется обычно в частотной области, реже во временной. Расчет нели- нейных видеосхем выполняется только во временной области, а расчет нелинейных радиочастотных схем — только в частотной. 1.3. Параметры и переменные электронных схем При проектировании схем приходится оперировать с различными параметрами и переменными. Все пара- метры можно разбить на две группы—параметры, схе- мы и параметры технического задания на проектирова- ние схемы. Параметры схемы разделяются на внутренние и внешние. К внутренним относятся параметры элементов (компонентов схемы) — транзисторов, диодов, микро- схем и т. д. Они подразделяются на электрофизические, топологические, электрические и режимные. Например, для МДП-транзистора концентрация примеси в подлож- ке — электрофизический параметр, длина и ширина ка- нала—топологические параметры, максимально допу- стимые напряжения и токи—режимные параметры. При объединении элементов в составе схемы все па- 16 . раметры элементов оказываются внутренними для схе- мы, а у схемы появляются новые внешние параметры. К внешним относятся быстродействие, помехоустойчи- вость, потребляемая мощность схемы и т. д. Внешние параметры, используемые для оценки качества работы схемы, называются выходными. Помимо выходных параметров-чисел схема и ее эле- менты имеют выходные характеристики, представляю- щие собой функциональные зависимости—частотные, переходные, амплитудные характеристики и т. д. Иногда выходные параметры и характеристики элементов и схем называют качественными показателями. К параметрам технического задания относятся пара- метры внешней среды (температура, давление, влаж- ность и т. д.), предельные режимные параметры (напря- жение питания и допуски на него, экстремальные зна- чения токов и напряжений в отдельных точках схемы), а также предельные значения выходных параметров и характеристик схемы. Техническое задание практически всегда выглядит как система одно- или двусторонних ограничений: г/^г/доп, у^г/доп, t/доп пип^г/г^г/доп шах, где Удоп—параметры технического задания, а у—парамет- ры проектируемой схемы. При проектировании схем используют также понятия управляемых внутренних параметров, изменение кото- рых физически легко осуществимо, и неуправляемых. В задачах оптимизации схем применяют понятие варьи- руемых параметров, изменяя которые можно улучшить значения оптимизируемых выходных параметров. Таким образом, варьируемыми параметрами могут быть толь- ко управляемые параметры. Кроме перечисленных параметров при проектирова- нии схем на ЭВМ используют понятие базисные пере- менные, т. е. переменные, относительно которых строит- ся система уравнений схемы (например, токи элемен- тов, узловые потенциалы и т. д.). Различают следующие основные типы базисов: пол- ный гибридный базис (образован токами и напряжения- ми всех элементов); полный однородный базис (обра- зован только токами или только напряжениями элемен- тов) ; сокращенный однородный базис (образован узловыми потенциалами или контурными токами); сокращенный гибридный базис (образован частью то- ков и частью напряжений элементов). Возможны и дру- гие базисы. 2—3469 17 1.4. Основные задачи проектирования При проектировании электронных схем приходится решать самые разнообразные задачи. Однако все зада- чи можно отнести к следующим четырем типам: расчет, анализ, параметрическая или структурная оптимизация, параметрический или структурный синтез. Расчет схемы — определение ее параметров и харак- теристик при неизменных значениях внутренних пара- метров схемы и постоянной структуре. Примерами за- дач расчета являются расчет статического режима, пе- реходного процесса, амплитудно- и фазочастотных характеристик. Анализ схемы—определение изменения выходных и режимных параметров схемы в зависимости от изме- нения ее варьируемых внутренних параметров. К зада- чам анализа относятся статистический анализ, анализ чувствительности, допусковый, температурный и другие виды анализа. Следует отметить, что граница между задачами расчета и анализа весьма условна, поэтому, например, статистический анализ и анализ чувствитель- ности часто относят к задачам расчета, а все задачи расчета называют одновариантным анализом. Отметим также, что процесс решения уравнений схемы часто на- зывают моделированием схемы. Параметрическая оптимизация — определение на- илучших в том или ином смысле значений заданных выходных параметров путем целенаправленного изме- нения варьируемых внутренних параметров (способ их изменения определяется конкретным алгоритмом опти- мизации). Структурная оптимизация — определение наилучших значений заданных выходных параметров путем целе- направленного изменения первоначально заданной структуры схемы, т. е. связей между ее элементами. Примером структурной оптимизации может служить введение дополнительных положительных или отрица- тельных обратных связей. Наиболее сложными являются задачи параметриче- ского и структурного синтеза. Параметрический синтез заключается в том, что для заданной структуры схемы определяются параметры ее элементов, обеспечивающие правильное функционирова- ние схемы, например ее работу в режиме триггера, а не двухкаскадного усилителя. В отличие от параметряче- 18 ской оптимизации при параметрическом синтезе пара- метры элементов схемы определяются не в процессе последовательного изменения некоторых первоначаль- ных значений, а путем расчета соотношений, в которых уже учитываются способ функционирования схемы и ее условия работоспособности. Схема, полученная в ре- зультате параметрического синтеза, не обязательно оп- тимальна, но обязательно работоспособна. Для улучше- ния ее выходных параметров можно использовать па- раметрическую оптимизацию. Если параметрический синтез позволяет сразу со- здать схему не только работоспособную, но и оптималь- ную, то такой синтез называется оптимальным парамет- рическим синтезом. В общем же случае параметриче- ский синтез можно рассматривать как процедуру определения начальной точки для последующей пара- метрической оптимизации. Все сказанное в отношении параметрических синте- за и оптимизации справедливо для структурных синтеза и оптимизации. Структурный синтез также можно счи- тать процедурой определения исходной структуры для последующей структурной оптимизации. 1.5. Виды расчета и анализа параметров схем Расчет схемы является простейшей задачей при ее проектировании. Расчет статического режима выполняется для по- строения карты режимов по постоянному току, а также определения различных статических параметров схем. Однако чаще всего статический режим рассчитывается для определения начальной рабочей точки перед после- дующим расчетом переходных процессов и частотных характеристик. Расчет переходных процессов служит основой опре- деления различных выходных динамических парамет- ров. Кроме того, он позволяет по виду переходного про- цесса в различных точках схемы составить заключение о работоспособности схемы. При расчете переходных и частотных характеристик проектировщик схемы должен указать интервал времени или диапазон частот, на ко- торых выполняется расчет. Расчет выходных параметров схемы основан на пред- шествующем расчете статического режима, переходных процессов или частотных характеристик. Обычно типы 2* 19 рассчитываемых выходных параметров не могут выби- раться проектировщиком схемы произвольно и должны соответствовать перечню выходных параметров, расчет которых может выполняться используемым комплексом программ. Чтобы снять это ограничение, в наиболее со- вершенных программах предусмотрена возможность за- писи любых выходных параметров по усмотрению про- ектировщика схемы с помощью средств входного языка (см. гл. 4). Следующей по сложности задачей проектирования является анализ схемы. В комплексах программ АСхП обычно предусматривается возможность выполнения следующих типовых видов анализа. Анализ параметрической чувствительности позволя- ет определить степень влияния изменения внутренних параметров схемы х на ее выходные параметры г. По- мимо самостоятельного значения эта задача играет вспомогательную роль в методах оптимизации, исполь- зующих производные. Статистический анализ, обычно выполняемый мето- дом Монте-Карло, позволяет найти практически любые статистические характеристики схемы. Результаты ста- тистического анализа часто представляются в виде ги- стограмм выходных параметров, по которым можно за- дать целесообразные границы отбраковки или, если эти границы уже заданы в техническом задании, оценить их жесткость. Кроме того, важным результатом статисти- ческого анализа служит оценка -вероятности годности схемы, определяемая как процент выхода годных схем, т. е. схем, по всем выходным параметрам удовлетворяю- щих техническому заданию. Для выполнения статистического анализа проектиров- щик схемы должен указать в задании на проектирова- ние, какие внутренние параметры варьируются, тип ста- тистического закона и значения его параметров для каждого варьируемого параметра схемы. При необхо- димости в ЭВМ вводится также корреляционная матри- ца варьируемых параметров. Анализ на наихудший случай является одним из ви- дов допускового анализа, при котором выходные пара- метры схемы рассчитываются для наихудшего сочета- ния внутренних параметров, определяемого в программе автоматически или задаваемого проектировщиком схе- мы на входном языке. Анализ влияния внешней среды состоит в расчете 20 зависимости выходных параметров г=г[х (/•)], где г— параметр внешней среды (температура, радиация и т. д.). Анализ выполняется автоматически в два эта- па: сначала рассчитываются зависимости x(r), a затем z{x). Типы зависимостей х(г) указываются проектиров- щиком в задании на анализ. Многовариантный анализ состоит в том, что проек- тировщик указывает несколько различных комбинаций (вариантов) параметров элементов схемы, а комплекс программ рассчитывает выходные параметры схемы, со- ответствующие каждому варианту, причем весь расчет выполняется за один вычислительный сеанс. В резуль- тате проектировщик получает результаты расчета сразу нескольких схем, отличающихся внутренними парамет- рами, и выбирает наилучший вариант. Режим много- вариантного анализа широко используется при проек- тировании схем на ЭВМ, во многих случаях заменяя режим оптимизации. Кроме того, задачи анализа на наихудший случай, статистического анализа, анализа влияния внешней среды также могут быть сформулиро- ваны как задачи многовариантного анализа. Для проектировщика представляет интерес вопрос о затратах машинного времени на выполнение различ- ных задач расчета и анализа. Однозначно ответить на него нельзя, так как время расчета схемы зависит, во- первых, от размера схемы и особенностей ее математи- ческого описания, во-вторых, от качества программы расчета, а в-третьих, от быстродействия ЭВМ, на кото- рой выполняется расчет. Экспертные оценки времени расчета приведены в гл. 5 (табл. 5.1). Время анализа схемы, очевидно, зависит от коли- чества расчетов, используемых в данном виде анализа. 1.6. Параметрическая оптимизация схем Оптимизация схемы — завершающий этап ее проек- тирования. При решении задач оптимизации с помощью программы АСхП главную роль играют три фактора— выбор критерия оптимальности, выбор алгоритма оп- тимизации и задание ограничений на внутренние, ре- жимные и выходные параметры схемы. Обычно в программе имеются библиотеки крите- риев оптимальности и алгоритмов оптимизации, из ко- торых проектировщик схемы выбирает наиболее под- ходящие по его мнению, для оптимизации данной 21 схемы. Ограничения на внутренние параметры схе- мы определяются проектировщиком исходя из усло- вий их физической реализуемости. Ограничения на ре- жимные и выходные параметры либо задаются проекти- ровщику заказчиком в виде технического задания, либо формируются самим проектировщиком на основе собст- венных представлений об особенностях работы и воз- можностях схемы. Решение задачи оптимизации схемы на ЭВМ закан- чивается после достижения выбранным критерием оп- тимальности своего экстремального значения. Очевидно, схема, оптимальная в смысле выбранного критерия, мо- жет быть неоптимальной в смысле других критериев. Существующие программы оптимизации допускают варьирование несколькими десятками параметров, что обычно достаточно для оптимизации схем. Время опти- мизации зависит от выбора начальной точки оптимиза- ции, критерия оптимальности, алгоритма оптимизации, времени однократного расчета схемы. В процессе опти- мизации расчет схемы выполняется, как правило, не менее нескольких десятков раз для любого алгоритма оптимизации, а допустимое общее время оптимизации ограничивается обычно несколькими часами. Поэтому практически на ЭВМ удается пока оптимизировать схе- мы, расчет которых требует не более нескольких минут, т. е. содержащие не более нескольких десятков транзи- сторов. Если размер схемы превышает допустимый для программы, то ее расчет и оптимизация выполняются по частям. Вообще, процессы получения оптимального варианта схемы отличаются большим разнообразием, в том числе активным участием проектировщика, и процедура опти- мизации, выполняемая на ЭВМ, часто служит лишь вспомогательным, а не окончательным средством опти- мизации схемы. 1.7. Структурный синтез схем Структурный синтез электронных схем — наименее формализуемая, и поэтому наиболее сложная для по- становки на ЭВМ задача проектирования. В отличие от оптимизации, для которой разработано огромное коли- чество алгоритмов, строгие алгоритмы синтеза структу- ры схемы весьма немногочисленны и разработаны глав- ным образом для линейных схем [4]. Поэтому главная 22 роль в структурном синтезе принадлежит человеку; его знания, опыт и интуиция определяют построение основ- ных связей и тем самым принцип действия схемы. По этой же причине в процессе синтеза широко использу- ются эвристические приемы, не имеющие строгого мате- матического обоснования, но позволяющие получить близкий к оптимальному результат, причем гораздо бы- стрее, чем строгими методами. Конкретные способы син- теза определяются типом схемы и ее размером. Типовая процедура формализованного синтеза со- стоит в том, что с помощью ЭВМ генерируются различ- ные варианты структуры, а проектировщик выбирает из них наиболее подходящие с его точки зрения. Спо- собы генерирования вариантов различны, но их можно разделить на три группы. В первой группе ЭВМ выдает проектировщику варианты известных структур, храня- щихся в ее памяти и, таким образом, играет роль элек- тронного справочника. Этот способ синтеза можно на- звать методом перебора известных решений. Во второй группе способов ЭВМ выдает модифици- рованную структуру, в основе которой лежит известная структура. Модификация может заключаться в том, что, например, в известную структуру вводятся допол- нительные отрицательные или положительные обратные связи, или в том, что в известной структуре схемы за- меняются ее отдельные элементы (блоки, каскады) раз- личными известными фрагментами. При этом сохраня- ются макроструктура схемы, например последователь- ность соединения фрагментов, но видоизменяется ми- кроструктура фрагментов. В памяти ЭВМ в этом случае должны храниться разнообразные варианты различных фрагментов схем. Наконец, в третьей группе ЭВМ генерирует произ- вольные структуры в том числе неизвестные ранее. Проектировщик сам или с помощью ЭВМ отбирает эти новые структуры. Этот способ синтеза называется син- тезом патентоспособных решений [5]. Помимо рассмотренных эвристических существуют строгие методы схемотехнического синтеза. Так, для цифровых схем по заданным таблицам истинности мож- но синтезировать сначала логические уравнения, а за- тем электрическую схему. В отличие от традиционного логического синтеза устройств на типовых логических схемах И—НЕ, ИЛИ—НЕ, в схемотехническом синтезе роль простейших логических элементов играют отдель- 23 ные транзисторы, сопротивления и т. д., что позволяет по заданному логическому выражению построить элек- трическую схему [6, 7]. ГЛАВА 2 ОСНОВНЫЕ МЕТОДЫ И АЛГОРИТМЫ АСхП В данной главе описаны основные алгоритмы, ис- пользуемые в программах АСхП при решении задач расчета, анализа и оптимизации схем. В связи с много- образием этих алгоритмов большое внимание уделяется их классификации. Ознакомление с ними позволяет обо- снованно выбирать состав алгоритмов при разработке программ, помогает пользователю квалифицированно задавать параметры алгоритмов при составлении зада- ний на расчет схем, способствует пониманию причин возможных остановок вычислительного процесса при расчете схем на ЭВМ. 2.1. Способы представления математических моделей схем В общем случае математической моделью схемы (ММС) называется система уравнений, формул, правил или любых других математически строгих соотношений, позволяющих найти с требуемой точностью определен- ные параметры или характеристики схемы. Понятие ММС трактуется весьма широко и привело к большому разнообразию способов их представления, отражающих различные степени детализации ММС. Приведем классификацию уровней представления ММС и соответствующие примеры. Уровень представления Пример ММС 1. Исходная модель- F(x)==0. (2.1) уравнение 2. Теоретическая мо- ХР+^ХР—^^ХР)}-^ дель-алгоритм (теоре- ><Р(л:Р). (2.2) тический метод) 3. Машинная модель- Схема на рис. 2.1. алгоритм 24 4. Структурная модель Схемы в гл. 7. программы 5. Языковая программная Текст программы на ал- модель горитмическом языке Назначение исходной модели—формализовать зада- чу, т. е. определить класс описывающих ее уравнений и принципиальную возможность их решения с помощью какого-либо теоретического метода или алгоритма (в данном случае—метода Ньютона). Теоретический ал- горитм реализуется далее с помощью машинного алго- ритма, назначение которого—учесть особенности чис- ленной реализации теоретического алгоритма. Следует отметить принципиальное отличие машин- ных алгоритмов от предполагающих идеальную точ- ность, устойчивость и сходимость вычислений теорети- ческих, заключающееся в том, что машинный 'алгоритм Начало Ввод начального приближения х^ и критерия точности вычислений е Вычисление F^x"), i- f,n\^——|e^'-a^r*—i Вычисление F'i(xf-^-AXj), i-1,n,j-f,n -, t=f,n, j~l,n , _, / F^^&X)-?^'')^'1 (Прощение матрицы Якоби (f(X)) •» 1————,————-j Вычисление F'(xp)~ff'(a!l'} ^ вычисление xr'^=я''-(F'(a!p))~i,j>i), (2.17) /К ik где элементы ;„• и ил берутся из крайнего правого столбца получившейся к i-му шагу части матрицы L и крайней нижней строки получившейся части матри- цы U. Таким образом, разложение по методу Краута — Холецкого соответствует схеме A-^i-^1)^-^-^2)- ..., а разложение по методу Дулитла — схеме l^l->h->^W->U^l2-^^m-> . . • Обратный ход в обоих вариантах выполняется оди- наково. Сначала решается система Lq=B но формуле -С) (] > 0. (2.18) а затем система Ux=q х / " ^.^LC-'» _ \i \ ' /='+i за ч •. г • "Ч--1 и. II- (2.19) Существуют и другие формулы для описания прямого LU-раз- ложения. Например, для 2-го варианта i_i / 1-1 -. i "t/ = "i/ — 2.1 lik^hi (J Ss: 0; 1ц = а^—З 'i^ki / "ii (i < ') • fc=l (2.20) В этих формулах не учитывается предыстория разложения, поэтому суммирование начинается от 6=1, одновременно отпадает необхо- димость пересчета матрицы А на каждом шаге. Преимуществом разложения по (2.20) можно считать возмож- ность рациональной организации вычислений суммы с минимальной потерей точности из-за ошибок округления при суммировании от- дельных слагаемых; в то же время предыдущий способ с пересчетом матрицы А на каждом шаге позволяет организовать, как будет по- казано далее, эффективный процесс расчета больших схем по частям. 2.2.3. Особенности машинных алгоритмов. Специфи- ка решения системы (2.4) состоит в необходимости вы- сокой точности решения и одновременно сохранения разреженности матрицы А. Для различения машинных алгоритмов, реализующих тот или иной теоретический алгоритм решения системы (2.4), определим следующие понятия. Метод упаковки ненулевых элементов — метод хра- нения матрицы А в виде списков — одномерных масси- вов, содержащих информацию о значении ненулевых элементов (НЭ) матрицы А и их номерах строк и столбцов. Критерий упорядочения для разреженности — способ определения оптимального порядка обработки уравне- ний в системе (2.4), при котором число новых НЭ (ННЭ) будет минимальным. Критерий упорядочения для точности — способ оп- ределения оптимального порядка обработки уравнений в системе (2.4) с целью повышения точности решения. Поскольку критерии упорядочения для разреженно- сти и точности могут приводить к разным последова- тельностям обработки уравнений, введем понятие пра- вила разрешения конфликтов между ними. Наконец, каждый критерий может быть реализован различными способами, поэтому вцедем понятие спосо- ба реализации критерия (в литературе часто критерий и способ его реализации ошибочно отождествляются). Рассмотрим вкратце различные способы реализации указанных понятий. Методы упаковки [14—25] можно разделить на уни- версальные, пригодные для матриц любого вида, и спе- 3-3469 33 j Методы упаковки Универсальные ' Специальные С произвольным С фиксированным порядком записи порядком записи Для полностью симметричных матриц Для структурно- симметричных матриц йез разделителей Со связными списками С номарзми строк к разделителями Для ленточных матриц Для блочных матриц С постоянной шириной ленть! С номерами столбцов и разделителями С переменной шириной ленты Рис. 2.5. Классификация основных методов упаковки разреженных матриц циальные, учитывающие структуру заполнения матрицы ненулевыми элементами (рис. 2.5). Все универсальные методы основаны на замене дву- мерного массива—матрицы А—одним или нескольки- ми одномерными массивами — списками. По содержа- нию списки можно разделить на списки значений НЭ, списки позиционных указателей, содержащие номера строк (столбцов) НЭ, и списки разделителей, содержа- щие информацию, позволяющую выделять в списке зна- чений НЭ отдельные группы, находящиеся в одной стро- ке (столбце). В качестве разделителей можно исполь- зовать число НЭ в строке (столбце), порядковые номе- ра (при чтении матрицы А подряд по строкам или столбцам) первых или последних НЭ и т. д. Во всех методах упаковки имеется список значений НЭ и хотя бы один из списков позиционных указате- лей — строчный или столбцовой. Вместо второго списка указателей часто для экономии памяти и удобства поис- ка элементов отдельных строк (столбцов) применяют список разделителей. Все универсальные методы упаковки можно разделить на две группы — с произвольным порядком записи, допускающим введение ННЭ без изменения индексации старых НЭ, и с фиксированным порядком записи, когда введение ННЭ не предусматривается, так как требует сдвига последующих старых элементов. Примером уни- версального метода упаковки с произвольным порядком записи мо- жет служить запись НЭ матрицы А в виде трех списков — значений НЭ, номеров строк и номеров столбцов. 34 К этим спискам можно добавлять информацию о любом числе ННЭ без изменения предыдущей информации. Если заменить один из списков позиционных указателей списком разделителей, например список номеров столбцов НЭ списком порядковых номеров первых НЭ в каждой строке, определяемых относительно первого НЭ в матрице, то получим метод с фиксированным порядком записи, так как введение ННЭ в k-ю строку потребует изменения в списке указателей — увеличения на единицу всех порядковых номеров, отно- сящихся к строкам, начиная с (k-{-\.)-u. Специальные методы упаковки ориентированы на оп- ределенную структуру матрицы А и описаны в [28]. Хотя экономия памяти безусловно играет важную роль при списочной записи матрицы А, однако нужно учитывать, что замена матрицы списками, требующими дополнительных вычислений для определения номера строки или столбца НЭ, замедляет поиск элементов в списках по сравнению с их поиском в исходном дву- мерном массиве. В связи с этим при выборе методов упаковки в пер- вую очередь нужно учитывать не занимаемую память (в этом отношении большинство методов примерно рав- ноценно), а скорость поиска НЭ в списках, для повы- шения которой может оказаться целесообразным вво- дить вспомогательные списки, например записывать диагональные НЭ, положение которых определяется лишь одним номером строки, в отдельный массив, что- бы не тратить время на их поиск в общем списке по общему правилу. Критерии упорядочения для разреженности доста- точно подробно изложены в [18—24]. В связи с тем что определение наилучшего порядка обработки уравнений при разложении А требует перебора огромного количе- ства вариантов, все критерии являются лишь квазиоп- тимальными, поскольку основаны на просмотре ограни- ченного числа вариантов. Основные характеристики этих критериев: длина просмотра /, т. е. количество уравнений, ана- лизируемых для однократного применения критерия; глубина просмотра h, т. е. количество шагов разло- жения, необходимых для однократного применения кри- терия; ширина просмотра g, т. е. количество различных ана- лизируемых порядков разложения на h шагах, необхо- димое для однократного применения критерия; мощность критерия w, т. е. количество уравнений в последовательности разложения, определяемое при З* 35 однократном применении критерия; в большинстве практических случаев ю=1; 1==п—i (i—номер шага разложения); ft=l, g==l. Приведем примеры наиболее распространенных кри- териев упорядочения: Кр 1: К/, если ИНЭ.<ИНЭ„ где ИНЭ„ ИНЭ, — число исходных (до разложения) НЭ в строках с но- мерами i и /"; Кр 2:i15—20; 3) при увеличении п (п> 100) различие в эффектив- ности критериев сглаживается. Последнее объясняется необходимостью многократ- ного применения критерия в так называемых неопреде- ленных ситуациях, когда критерий удовлетворяется одновременно для разных номеров строк (пар «стол- бец—строка») и, следовательно, порядок разложения может выбираться неоднозначно. Для устранения неод- нозначности нужно увеличивать параметры hug кри- терия, но это ведет к дополнительным затратам време- ни и усложняет алгоритм реализации критерия. 36 Способы, реализации критериев упорядочения для разреженности. В связи с изложенным целесообразно использовать те критерии, реализация которых наибо- лее проста. На практике получили распространение два критерия: Кр 2 и Кр 4. Подсчет числа ННЭ в Кр 2 легко осуществляется на основе анализа формул разложения. Например, при прямом LU-разложении на каждом шаге разложения количество ННЭ в пересчитанной оставшейся части ма- трицы А зависит от того, какая пара «столбец — стро- ка» была выбрана на предыдущем шаге разложения. Очевидно, следует выбирать ту пару, при разложении по которой образуется минимальное число ННЭ. Из (2.17) видно, что при пересчете матрицы А на месте нулевого элемента а ^"^ появится ненулевой элемент а'0 если одновременно /д=/=0 и Uih.=^Q {j>i, k>i). Отсюда следует, что для минимизации числа ННЭ нужно выбирать на каждом шаге разложения ту пару «столбец — строка», в которой минимально число не- нулевых произведений /,„ иц,- Если после очередного шага разложения несколько пар «столбец—строка» да- ют одинаковое минимальное число ННЭ, то выбирается пара, содержащая максимальное число НЭ, так как при этом из дальнейшего рассмотрения исключается максимальное число элементов. Выбрав по этому пра- вилу последовательности рассмотрения пар «столбец— строка» матрицы А, нужно перенумеровать уравнения и сформировать матрицу А для новых номеров. Этот спо- соб обычно называется критерием Берри [21]. Приме- нение критерия Берри при решении, например, раз- реженной системы (2.4) из 100 уравнений с коэффици- ентом заполнения йз=4% ускоряет решение почти в 40 раз. В табл. 2.1 приведены сравнительные результаты расчета некоторых схем по программам, использующим данный критерий и без него. Как видно из таблицы, за- траты памяти и времени при использовании критерия упорядочения почти пропорциональны числу узлов п^ т. е. порядку решаемой системы уравнений. Для реализации Кр 4 удобно использовать способ Марковица [22], в соответствии с которым для оче- редного исключения выбирается та неизвестная, коор- динаты которой f, j удовлетворяют условию Л", 37 Таблица 2.1 LU -разложение без LU -разложение Число учета разреженности с учетом разреженности Тип схемы узлов в схеме Число Число Число Число операций ячеек операций ячеек Мультивибратор 11 553 121 106 51 10-звенная цепь 22 1507 484 208 120 Схема на МДП- 35 14245 1225 570 222 транзисторах Тест-схема 100 333 000 10000 2264 807 где ;Vi, N, — число НЭ в строке i и столбце /'. Величи- на Ni, равна максимально возможному числу ННЭ при дальнейшем разложении в самом неблагоприятном слу- чае, когда каждый ненулевой элемент столбца порожда- ет ННЭ во всех возможных позициях соответствующей строки. Этот способ упорядочения проще предыдущего, мало уступая ему по эффективности в смысле числа ННЭ [18]. Заканчивая рассмотрение критериев упорядочения и способов их реализации, можно отметить, что при всех различиях используе- мых критериев в результате их применения к большим системам уравнений исходное хаотичное расположение ненулевых элементов в матрице А обычно приобретает вид, изображенный на рис. 2.6 [3, 18, 26]. Очевидно, это не случайно. Объяснить это явление мож- но, опираясь на следующую теорему [9]. Если для некоторого /'(t) элементы матрицы А удовлетворяют условиям а,,=0, 1=1, ..., р< е может не выполняться сколь угодно долго. В связи с этим вычислительный процесс при реше- нии (2.4) нужно строить с учетом требований миними- зации ошибок на каждом шаге решения. Обычно для этого применяют какие-либо методы главных (ведущих) элементов. Главным называется элемент матрицы, на который выполняется деление на очередном шаге пря- мого хода решения системы (2.4). Деление на макси- мальный главный элемент дает минимальную ошибку округления, поэтому на каждом шаге прямого хода предпочтение отдается строке (паре «столбец — стро- ка») с максимальным значением диагонального эле- мента. Главный элемент может отыскиваться среди либо только диагональных элементов, либо всех элементов матрицы А (поиск по всему полю матрицы). В первом 39 случае нужно переставлять и строку и столбец, а во втором — только строку или столбец, чтобы выбранный элемент стал диагональным. Правило разрешения конфликтов между критерия- ми упорядочения для разреженности и для точности со- стоит в следующем. Устанавливается допустимое минимальное значение s диагонального элемента a,i(i==l, ..., п). Если все а,ц больше е, то работает критерий упорядочения для разреженности, если для некоторых а„ это правило не выполняется, то приоритет отдается критерию упоря- дочения по точности, т. е. строка с минимальным ац об- рабатывается при LLJ-разложении или гауссовом ис- ключении в последнюю очередь. Реализация этого пра- вила требует переупорядочения уравнений по различ- ным критериям. Чтобы избежать этого, иногда исполь- зуют метод диагональной модификации [3, 38]. 2.3. Решение систем нелинейных уравнений 2.3.1. Основные проблемы. Исходная модель-уравне- ние этого класса задач имеет вид F(x)==0, (2.23) где F(x) —вектор-функция, х—вектор. При решении (2.23) на ЭВМ возникают следующие основные проблемы. 1. Численные методы, гарантирующие решение (2.23), требуют выполнения определенных условий, на- пример квадратичности F(x), ограниченности нормы ilF^x)!! известным значением и т. д. При расчете элек- тронных схем эти условия обычно неизвестны и гаран- тии сходимости нет. Поэтому первой проблемой реше- ния (2.23) на ЭВМ является обеспечение высокой алго- ритмической надежности решения, т. е. сходимости в возможно более широкой области начальных прибли- жений, при возможно более разнообразном виде функ- ций F(x) ив возможно более широком диапазоне чис- ленных значений ее параметров. 2. Отдавая приоритет требованию алгоритмической надежности, необходимо решать (2.23) с минимальными вычислительными затратами. К сожалению, численных методов, одновременно удовлетворяющих требованиям высокой алгоритмической надежности и малых вычис- лительных затрат, пока не существует. Поэтому второй 40 /•'(-V) =0 Прямые методы Методы оптимизации Метод простых итераций Методы нулевого порядка Метод Ньютона Методы первого поря Метод Ньютона — Бройдена Методы второго порядка Методы продолжения решения по параметру С ограниченным параметром С неограниченным параметром (методы установления) С одной итерацией Явные Неявные Метод движущейся области сходимости Рис. 2.7. Классификация основных методов решения систем нелиней- ных уравнений F(x)=0 основной проблемой решения (2.23) является обеспече- ние компромисса между алгоритмической надежностью и вычислительными затратами. 2.3.2. Основные теоретические методы. Основные ме- тоды решения (2.23) можно разделить на две группы (рис. 2.7): прямые методы и методы оптимизации. Среди прямых методов наиболее часто используется метод Ньютона (2.3). Его достоинства состоят в том, что при решении систем (2.23), описывающих реальные схемы, метод Ньютона, во-первых, всегда имеет область сходимости и, во-вторых, обеспечивает высокую ско- рость сходимости (квадратичную по сравнению с ли- нейной для метода простых итераций). Недостатками метода Ньютона являются малый размер области сходи- мости (обычно меньше, чем для метода простых ите- раций) и необходимость задания точки начального при- ближения достаточно близко к точке решения — в про- тивном случае скорость сходимости заметно снижается. Метод Ньютона эффективен, если F(x)—выпуклая функция, в этом случае решение с ошибкой е== ^lO^—'lO-4 достигается за 3—5 итераций. 41 Теоретические условия сходимости метода Ньюто- на—существование и ограниченность обратной матрицы Якоби в точке начального приближения, ограничен- ность и непрерывность первых и вторых производных •функций [30]—носят характер теоремы существования и не могут быть использованы на практике для опреде- ления положения и размера области сходимости. В связи с этим начальные условия х101 в общем слу- чае приходится задавать эмпирически; при этом сходи- мость уже не гарантируется. Для повышения надежно- сти сходимости можно использовать модификацию ме- тода Ньютона—метод Ньютона—Бройдена [31] F/(xP)AxP=—rpF(xP), AxP^xP+'—xP, (2.24) где Гр — параметр, выбираемый на каждой итерации из условия минимума (р+1)-й невязки: F(xP+l)= =F[x?—rp(F/(xP))-lF(xP)]—«-min. Минимизация невязки г Р(х?+1) на каждой итерации расширяет область сходи- мости этого метода по сравнению с методом Ньютона. Методы оптимизации основаны на эквивалентности результатов решения системы (2.23) прямыми методами и минимизации какого- либо из функционалов •Ф,(х)=3р2, (х); Ф, (х) =3|F, (x)|; Фз(х)=тах|Рг(х)| i i ' (2.25) хорошо отработанными методами оптимизации. григг [29]. При этом количество решений системы (2.23) равно количеству минимумов в (2.25). Отсюда следует, что локальные методы опти- мизации можно использовать для расчета статического режима схем с одним состоянием равновесия (ключи, усилители и т. д.), если же схема имеет несколько состояний равновесия (триггеры), то локаль- ные методы могут .и не дать нужного решения |29J. В целом методы оптимизации не нашли широкого применения для решения системы (2.23) при расчете схем, так как в большин- стве своем они уступают методу Ньютона по скорости сходимости, не обеспечивая в то же время высокой надежности сходимости. Повышение скорости сходимости этих методов до уровня метода Ньютона требует вычисления матрицы Гессе Ф"(х), что неприемле- мо при расчете схем на ЭВМ из-за больших затрат времени. 2.3.3. Повышение надежности сходимости. Для по- вышения надежности сходимости методов решения си- стемы (2.23) обычно используют принцип продолжения решения по параметру, служащий основой построения большинства методов, классификация которых приведе- на на рис. 2.7. Рассмотрим группу методов с ограничен- ным параметром. В этих методах система F(x)==0 за- меняется эквивалентной системой F(x, г)=0, где пара- 42 метр re. (/"min, Гюах) выбирается так, что решение Хо си- стемы F(x, rmin)==0 известно, a F(x, Гщах) ==F(x). Для удобства изложения введем оператор algF(x), обозначающий результат одной итерации решения си- стемы F(x)===0 любым численным методом, и оператор AlgF(x), обозначающий результат окончательного ре- шения той же системы [29]. Тогда метод продолжения решения по параметру с ограниченным параметром и одной итерацией на каждом q-u продолжении (шаге). кроме последнего, можно записать в виде x„+i==AlgF(x„, /-,,+,), <7==0, 1, ..., п - 1; (2.26а) r^=W; r,=^,„; г„=г^; (2.266) x*=x„=AlgF(x„-„ r^j, (2.26в) X где х* — решение системы; q— номер продолжения. Таким образом, в этом методе однократное решение системы (2.23) заменяется серией итераций с изменяю- щимся параметром на каждой итерации, а на последнем продолжении вместо одной выполняется полный цикл итераций. От этого метода отличается метод движущейся об- ласти сходимости [29] x,+,==AlgF(x г +,), q==0. I, ..., й- 1; (2.27а) Y (2.276) В этом методе, если Хд^бХд, где бХд — область схо- димости около точки \q, то возмущение параметра Лг,,= ==rg+i—Гд передвигает область сходимости и выполня- ется x^+is6xg+i (рис. 2.8). Условием надежности мето- да является XgSfiXg+i, так как в этом случае при выбо- ре в (2.27а) \q в качестве начального приближения ре- шение (2.27а) не расходится. Это условие легко выпол- нить выбором Гд+i. Конкретные машинные алгоритмы реализации методов (2.26), (2.27) отличаются типом и /^^\ /'"^'^ ^^\ i (• ).'...' i .} .\ ... i' (.!.' s^y v^w V^y ^-\ iFSn-i /fxn-e'x' /'"^!-/ '~Л~ГуЛ-С '(.х^тах)" F(S) -\ ifx^i r-r^ Рис. 2.8. Иллюстрация метода движущейся области сходимости 43 способом изменения параметра г и функции rq+\=f(rq). Обычно в качестве г выбирают напряжение питания Е, причем Eq+i==Eq+^E, где A? подбирается на каждом •продолжении из условия сходимости в (2.27а). Метод (2.26) более быстрый, чем (2.27), но менее надежный, так как Лгд выбирается произвольно (неза- висимо от Sxq+i), и к моменту решения (2.26в) может оказаться, что Xn-i^6x„. Рассмотрим теперь группу методов с неограниченным изменением параметра гг(0, оо). В этом случае вме- сто (2.23) можно рассматривать систему F(dx/dr, x) == ==0, в частности систему dx(r)/dr=F(x(r)) (2.28) при условии, что dx/dr-fO при v->-oo. Если (2.28) составляется не для всех / уравнений, входящих в (2.23), а только для некоторых, то можно записать систему дифференциальных и конечных урав- нений dxi(r)ldr=Fi(x(r)), t=l,...,^, (2.29a) F,(x)=0, j==k+\, .... I, (2.296) которая решается совместно. При этом систему (2.296) можно решать одним из предыдущих способов. Пусть k=l, тогда, решая (2.28) явным методом Эйлера, получаем x„+i=Xn+ArF(Xn). Легко видеть, что этот алгоритм соответствует решению (2.23) методом простых итераций. Тогда условие устой- чивости явного метода Эйлера Лг^1/|^[ шах, где |^]max— максимальное собственное значение матрицы F(x), можно рассматривать как условие сходимости метода простых итераций. Известно, что для плохо обусловленных матриц F^x) явный метод Эйлера из-за ограничений на &.г сходится очень медленно. Для ускорения решения целесообразно использовать неявный метод Эйлера, не имеющий ограничений на шаг, x^i=x„+^F(x?„+i). Здесь при каждом п производится вычисление Xn+i методом про- стых итераций. Для повышения надежности сходимости вместо простых ите- раций целесообразно использовать метод Ньютона. В этом случае нет необходимости получать явную систему (2.28), достаточно ре- шать неявную систему 'F(d\fdr, х(г))=-0, в которой заменить dx/dr разностным оператором. В простейшем случае для неявного 44 метода Эйлера нужна замена d\/dr== (Xn+i—х„)/Д/"; тогда алго- ритм будет иметь вид F^xPn+i, Xn, Ar)'AxPn+i=—F(xPn+i, x„, Дг). Здесь Дг можно выбирать достаточно большим. В конкретных машинных алгоритмах, реализующих методы решения систем (2.28), (2.29), в качестве пара- метра г обычно выбирается время t, а в качестве x(r)— реактивные переменные itc(t), t'i.(0, тогда условие dx(r) /dr\->-0 эквивалентно условию ic(t)->-0, иь(1)->0, что автоматически выполняется в статическом режиме. 2.3.4. Особенности машинных алгоритмов. При реше- нии системы (2.23) машинными алгоритмами основные особенности заключаются в способах задания началь- ного приближения х(°) и контроля сходимости. Выбор начального приближения определяется толь- ко требованием обеспечения сходимости. Значение х^ либо задается пользователем на основе знания им при- близительного решения системы (2.23), либо выбирает- ся в программе автоматически, как в методе движущей- ся области сходимости. Контроль сходимости при решении (2.23) методом Ньютона (2.3) выполняется одним из двух способов — по норме поправки llx^+^—x^ll^lAxfP)!! или по норме невязки ||F(xtP))]|, р — номер итерации. Если JiAxtP^I^ =0, то и ||F(x

)||==0, однако если ||F(x

) ||<е, то мо- жет быть, что ЦЛх^Ц^е, и наоборот. Поэтому целесо- образно в качестве критерия сходимости потребовать max(l|Ax(P)||, ||F(x tn—m) , (2.37) где бп(-), 6n+i(-)—дискретизированные (численные) производные в точках tn, tn+i- Выражая бп или 6n+i через Хп+\, Хп, •.., Хп-т и раз- решая получившееся уравнение относительно л-n+i, по- лучаем окончательную форму решения нормальных ОДУ — явным способом A'n+l==:•' (Хп, Xn—l, • • ., Хп—т, tn, Гта—1, • • •i tn—m) (Z.do) и неявным способом Xn+l=^г(Xn+i, Хп, . . ., Хп—т, tn+l, tn, . • ., tn—mi. (2.39) Для неявных форм (2.33), (2.35) различные способы дискретизации производных и интегралов приводят в конечном счете к уравнению вида Р (dn+lXn+l, ЧпХп, • • ., Чп-тХп-т, Pn+lSn+l^", '^n^nX, . . . . . ., рп,_,пбп-тп^> tn+l, tn, • • •> tn—m)==0. (2.40) При ап-ц=0, pn+i==0 это уравнение соответствует явным методам решения (2.35), а при ctn+i^O, pn+i^ ^О—неявным. Явные методы допускают обычно анали- тическое выражение Xn+i в виде (2.38), в неявных ме- тодах для определения Xn+i в общем случае нужно при каждом га+1 решать (2.40) каким-либо численным ме- тодом до сходимости, т. е. ^+i=AlgF(^„+i,^, ...,^„-^,^+1, ...,^-J. (2.41) •~П\-1 Заметим, что распространенная форма дискретиза- ции производной {)х== (Xn+i—Хп) I \t имеет одинаковый вид для n-v. и (п+1)-й точек, т. е. 6n+iX=6nX. Поэтому тип метода (явный или неявный) при использовании этой дискретизации определяется только значением ата+i. Так, в уравнении F(x-n+\, (Хп+\—Xn)/&t)=0 заложен неявный метод, а в уравнении F{Xn, (Xn+i—Хп) /ЛО = ==0 -- явный. Приведенные соображения позволяют для неявной формы ОДУ (2.40) определить явный или неявный тип заложенного в них метода решения уже на этапе со- ставления исходной модели-уравнения. 2.4.2. Основные характеристики численных методов решения ОДУ. К основным характеристикам численных 47 методов решения ОДУ относятся устойчивость и точ- ность. Устойчивость. При численном решении ОДУ на каж- дом шаге вместо точного решения x{tn) получается приближенное решение Хп. Пусть начальное условие Xo=x(to). Тогда на первом шаге получим, что Хг^- ^x^i), во-первых, из-за ошибок округления e°i при вычислениях и, во-вторых, из-за методической погреш- ности самого численного метода e"i, например погреш- ности дискретизации производной или замены функции f(x, t) полиномом Р(х, t}. На втором шаге получим, что Xt^=x(ti} не только из-за ошибок е°2 и s"2, возникаю- щих по указанным причинам, но и потому, что Xi== =(f(x\), a x\ было получено с ошибкой. Ошибка из-за использования при вычислении Хц+i неточных предыдущих значений Хп называется переход- ной ошибкой e^+i. Таким образом, если е\-=х\—л'(-оо. Уменьшение локальной ошибки ета требует соответственно выполнения условий е^—^-О при /г-»-0 (тре- бование аппроксимации, причем если е"^0 (h^, то го- ворят об аппроксимации порядка /Iй) и в"-»-0 при п-^оо (требование устойчивости). Таким образом, из аппрок- симации и устойчивости следует сходимость метода. Очевидно, если метод устойчив, то при п->-оо е" затуха- ет и локальная ошибка en ограничена и определяется только методической погрешностью, не зависящей от п. Если при п->оо ошибка еп-»-°о, то метод называет- ся неустойчивым, а если en растет, но не быстрее, чем 48 само решение Хп, то — относительно устойчивым (огра- ничена относительная погрешность вп/^та). Отметим, что, если даже методические ошибки в" отсутствуют, не- избежны ошибки округления е°, которые при отсутствии устойчивости приведут к расходимости метода. Для количественного определения условий, при которых метод устойчив или неустойчив, в теории численного решения ОДУ вво- дится понятие абсолютной устойчивости метода, согласно которому разностный метод называется абсолютно устойчивым, если все кор- ни г, соответствующего ему характеристического уравнения Р(г)== =0 находятся в круге единичного радиуса. Обычно в качестве тестового уравнения при анализе устойчи- вости метода используется линейное уравнение х=Ах или соот- ветствующее ему y='ky, где у^Мх, \=МАМ~^—диагональная- матрица, элементами ^.i^di-^jfui которой являются собственные значения матрицы Л. Тогда характеристическое уравнение можно записать в виде Р(г, 1Л) =0, а границу условия абсолютной устой- чивости рассматривать либо как окружность единичного радиуса в комплексной плоскости Rez, Imz, либо как соответствующую этой окружности кривую в комплексной плоскости Re/iA.=/io', Im/i/.== =jh(u. Уравнение этой кривой получается, если в характеристиче- ском уравнении Р{г, /Л)=0 положить |г|==1 и решить его относи- тельно /Л. Полученная кривая разделяет комплексную плоскость К на области абсолютной устойчивости и неустойчивости. Вид этих, областей порождает различные определения абсолютной устойчи- вости. Одним из наиболее общих является определение Л (а)-устойчи- вости (32, 33], согласно которому метод Л (а)-устойчив, если для системы линейных ОДУ с матрицей Л, решаемой численно с по- стоянным шагом h, все собственные значения ^^а'»-)-/®; матрицы hA лежат внутри угла S(d, —а), 0<а<я/2, принадлежащего ле- вой полуплоскости (рис. 2.9,а). Частными случаями Л (а)-устойчи- вости являются Л-устойчивость (во всей левой полуплоскости) (рис. 2.9,6), Л (0)-устойчивость (при а->-0), Ло-устойчивость (на отрицательной полуоси). Помимо Л (а)-устойчивости иногда исполь- зуют понятия жесткой устойчивости (рис. 2.9,в), при которой метод должен быть Л-устойчивым в области Ri и точным в области Rs, и ограниченной устойчивости (рис. 2.9,е), когда метод устойчив лишь в строго ограниченной области значений >.,h. Re/?A S) Рис. 2.9. Области устойчивости численных методов решения ОДУ: а—для Л (а)-устойчивого метода; б—для Л-устойчивого метода ((х==я/2); в—для жестко устойчивого метода: г—для ограничен- но устойчивого метода 4-3469 49 Необходимо подчеркнуть, что все приведенные опре- деления типов устойчивости справедливы только для линейных ОДУ, решаемых с постоянным шагом. Теория устойчивости решения нелинейных ОДУ с переменным шагом не разработана, имеются только частные резуль- таты, подтверждающие, например, Л-устойчивость не- явного метода Эйлера. Однако практика показала, что выводы, полученные на основе теории устойчивости ли- нейных ОДУ, можно использовать для нелинейных. Точность. Для характеристики точности численного метода решения ОДУ обычно применяют понятие «по- рядок метода», которое для разных групп методов име- ет разный смысл. Так, для методов, основанных на за- мене точного решения отрезками ряда Тейлора, порядок метода означает порядок производной или степень А" в первом отброшенном члене ряда. Порядок 0 (/г2) озна- чает, что первый отброшенный член равен x"(t)h2/2. В методах Рунге — Кутта обычно (но не всегда) поря- док метода равен числу этапов для выполнения одного шага. В линейных многошаговых методах (ЛММ), ос- нованных на замене точного решения интерполяционным полиномом, порядок метода равен степени интерполи- рующего полинома. Рекомендации по выбору численных методов реше- ния ОДУ сводятся к выбору метода, обеспечивающего максимальную скорость решения при достаточной точ- ности и безусловном выполнении требования устойчиво- сти. Следует отметить, что с повышением порядка ме- тода, т. е. точности, его устойчивость по сравнению с методом того же класса, но более низкого порядка ухудшается. Классификация численных методов реше- ния ОДУ приведена на рис. 2.10. В отношении устойчивости этих методов можно от- метить следующее: доказано (теорема Далквиста), что линейные неяв- ные методы выше второго порядка не могут быть Л- устойчивыми, но могут быть Л (а) -устойчивыми; с увеличением порядка Л (а)-устойчивого метода об- ласть его устойчивости уменьшается; так, для Л (а) - устойчивых методов типа ФДН (формул дифференциро- вания назад) значение а уменьшается от л/2 рад для методов 1-го и 2-го порядков до а==0,328 рад для ме- тода 6-го порядка [24]; не все неявные методы обладают свойством Л(а)- устойчивости, например его нет у неявных методов Рун- 50 те — Кутта выше второго порядка, неявных методов Адамса (методов Адамса—Мултона); все явные методы не имеют области Л (а)-устойчи- вости и относятся к ограниченно устойчивым, т. е. об- ласть их абсолютной устойчивости всегда ограничена (см. рис. 2.9,г). Для точек на границе области (на оси Re/Л) можно записать /irp|X,i|max==C, откуда следует ограничение на шаг интегрирования явными методами /г^гр=С/|?..|шах, (2.42) где С—постоянная данного метода. Так, для явного метода Эйлера (Xn+i=Xn-\-hf(Xn, tn)) С=\, для явного метода Рунге—Кутта 4-го порядка С =2,78. Если схема описывается линейным ОДУ вида х==Ах, то действительные собственные числа Ki матрицы Л связаны с постоянными времени т, схемы соотношением Ti==l/^i. Отсюда с учетом (2.42) следует, что для так называемых жест- ких схем, т. е. схем с большим разбросом постоянных времени innaxS>Timin, метод решения ОДУ должен быть устойчив в диапа- зоне шагов /1щ1п3 в силу рассмотренных выше причин малонадежна из-за большой дисперсии. Контроль потери устойчивости решения ОДУ осу- ществляется методами, аналогичными или близкими ме- тодам контроля ошибки. Очевидно, поскольку рост ло- кальной ошибки е предшествует потере устойчивости, то контроль ошибки предотвращает потерю устойчиво- сти. В неявных Л-устойчивых методах под потерей устой- чивости понимается расходимость итераций при расчете Хп+\- Формальным признаком потери устойчивости слу- жит несходимость итерационного процесса за заданное число итераций с заданной ошибкой е. В явных методах потеря устойчивости означает пре- вышение допустимого шага. Формальным признаком по- тери устойчивости обычно служит появление периодиче- ских осцилляции («гребенки») значений Хп, Xn+i и т. д. В обоих случаях для сохранения устойчивости нуж- но уменьшать шаг h. В неявных методах при этом вновь начинает соблюдаться условие ^(0) ебл:п+1, а в явных— условие /K/irp. Выбор шага осуществляется на основе результатов контроля точности или устойчивости. В неявных методах чаще всего используется следую- щий алгоритм изменения шага. Назначаются оценки Бтш И бтах ДЛЯ ЛОКаЛЬНОЙ Ошибки 6. Тогда 56 /1нов==/1ст при ещ1п<е<етах; /1нов>Йст ПрИ е-^emln", йновешах. (2.47) Конкретное соотношение между йнов и her зависит от способа контроля локальной ошибки. Если контроль выполняется прямым способом, с оценкой значения ошибки 8k, n+i, например, по (2.46), то новый шаг вы- числяется по формуле йков=^/е,.^Л„ (2.48) где едоп==ешш, если e,k,n+i < emin; едоп=вшах, еслие„,п+1> >Бтах. Если контроль локальной ошибки выполняется косвенным способом, например по числу итераций р, требуемых для вычисления Хп+\ с ошибкой не более вдоп, то вместо (2.47) можно использовать правило h»oB==hcT При pmW^P^Pmnx; hy.o-&=Ch^ при pJOmax. (2.49) Здесь С — коэффициент изменения шага, субъектив- но задаваемый пользователем, обычно С=2. Значения ртт, ртах зависят от алгоритма расчета Xn+i; обычно д^.я алгоритма Ньютона выбирают pmin=2—-3, ршах== =-5—6. В явных методах алгоритм определения шага фор- мулируется следующим образом: /гнов=^ах "Р11 Д-^Д-^оп; ДХд (2.50) "• /1ст при Дх > Д-^оп. где Ал'доп — задаваемое предельное приращение л: на один Шаг, ^Х=Хп—Xn-l, hcT==tn—tn-l, /1нов'=^эт+1—tn', "max— предельно допустимый шаг, задаваемый пользователем, Ограничение на максимальный шаг hma^. задается пользователем, исходя из особенностей решения кон- кретной задачи так, чтобы не пропустить быстрых из- менений решения x(t), например при определении фрон- тов импульсов. Иногда в неявных методах вводят ограничение на минимальный шаг, чтобы в случае невозможности обес- печить заданную локальную точность, ограничить про- цесс уменьшения шага и при ft) входит как слагаемое собст- венная проводимость у„ м-полюсника, если f-й узел схе- мы одновременно является полюсом га-полюсника. Рассмотрим пример. Любой биполярный транзистор можно рас- сматривать как трехполюсник с матрицей проводимостей э Э Уээ Утр= б —Уэб к —Уж —Уэб —Уэч Убб —'/кб —Укб Укк (2.61) Формулы для элементов этой матрицы t/ээ, Уэб и т. д. различны и зависят от выбранной схемы замещения транзистора (см. гл. 3). Рассматривая транзистор в схеме иа рис. 2.14 как п-полюсник, мож- но сформировать матрицу узловых проводимостей всей схемы Y((p) 1 2 34 1 ^ 1 ~^ + Ул -У я. о о 1 1 2 -Уа Уд+ ~R, +кГ +У06 -Укб —Уэб 1 3 о -Укб R. •+ Укк —Уж 4 о -Уэб —Ум Уээ+-^ Подводя итог рассмотренным методам формирова- ния ММС в базисе узловых потенциалов, можно сфор- мулировать следующий алгоритм формирования ММС. 1. Уравнение каждого элемента схемы представляет- ся в виде 1э==Уэ(р+1э, где Уэ, 1э—матрица полюсных проводимостей элемента и вектор полюсных токов. 2. Выделяются и обнуляются двумерные или одно- мерные массивы для записи в них матрицы узловых проводимостей схемы Y и вектора узловых токов I, яв- ляющихся основными компонентами ММС. .64 3. Рассматривается очередной k-и элемент схемы. Компоненты его вектора 1э и матрицы \э заносятся на соответствующие позиции массивов I, Y и суммируются с их содержимым. 4. Если список элементов не исчерпан, то повторя- ется п. 3, если исчерпан, то формирование ММС на этом заканчивается. 2.6.3. Учет уравнений реактивных элементов. Рас- смотренные выше случаи относились к ММС для рас- чета статического режима. Если формируется ММС для расчета переходных процессов, то в составе вектора 1((р) и матрицы Y((p) нужно учесть уравнения (2.31)— (2.32) реактивных компонентов—емкости и индуктив- ности. Для этого выбирают некоторый интервал А^ и эти уравнения дискретизируют на интервале М по форму- лам численного дифференцирования или интегрирова- ния. Поскольку в методе узловых потенциалов уравне- ния элементов должны иметь вид i==f(u), то дискрети- зированные уравнения реактивных элементов представ- ляются в виде ic=ycUc+Ic; lL=УLU•L-\гIL, (2.62) где ус, UL—проводимость реактивных ветвей; Ic, IL— источники тока. Конкретная формула численного дифференцирова- ния или интегрирования определяет тип метода решения ОДУ—Эйлера, трапеций, ФДН к т. д., а точка, в кото- рой выполняется дискретизация, — класс метода — яв- ный или неявный. Пусть уравнения электрического равновесия (балан- са токов и напряжений) составляются для точки tn. Тогда все формулы дискретизации производных ;'с== ==Cdiic/dt; iiL==Ldii/dt, приводящие к выражениям вида ^C,n+[=УcUc,1l+l-\-Ic(tn, . . ., fn+1-m), lL,n+l==l/LUL,n+l+lL(tn, . . ., tn+l-m), (2.63) соответствуют неявным методам решения ОДУ, а при- водящие к выражениям вида lC,n=УcUc,n+\-{-Ic(tn, ..., tn+\-m), lL,n+l==ybUL,n-{-lL(fn, ..., tn+l-m) (2.64) —явным методам. Значения токов источников /с, IL 5-3469 65 определяются через значения ис, ic и UL, IL в моменты tn, • • -, tn+1—ffl. Рассмотрим пример. Простейшая дискретизация производных t'c.n+l^C^c.n+l——Uc,n) /At', UL;n+l'wL(iL;n+l——lL,n)/^t соответствует неявному методу Эйлера. Приведем ее к виду (2.63): ?- ^с.п ^C.n-^•\ = Ы. UC, п+1 — Д^ ' д< , . ^•L.n^•\ =•~[~UL.^t+\ + t•L,n\ отсюда видно, что yc=CfAt; yL=At/L; lc=Cuc,n/^t; /L==i'i,,n. При формировании ММС по методу узловых потенциалов про- водимости ус, уь независимо от метода решения ОДУ (явного или неявного) включаются по описанным выше для многополюсников правилам в матрицу узловых проводимостей V. Если метод решения ОДУ — неявный, то в вектор узловых токов I (<р) заносится вели- чина i'c,n+i. или ii,,n+i> вычисленная по (2.63) для момента tn-f-i. Если метод решения ОДУ — явный, то в вектор 1 (ср) заносятся толь- ко источники /с, IL, вычисленные по (2.64) для момента t-n. Более подробно реализация явного метода узловых потенциалов изложена далее в п. 2.8.2. « 2.7. Алгоритмы расчета статического режима В основе алгоритмов расчета статического режима лежит решение уравнения (2.57). При решении (2.57) чаще всего используется алгоритм Ньютона (2.58), од- нако он малонадежен. Остановка расчета при реализа- ции алгоритма (2.58) на ЭВМ происходит вследствие переполнения разрядной сетки, возникающего либо из- за расходимости вычислений, либо даже в сходящемся в принципе процессе из-за появления слишком больших чисел, например при вычислении экспонент в уравнени- ях токов р-п переходов (рис. 2.15,а, б). Ниже приво- дятся алгоритмы решения (2.57), отличающиеся высо- кой надежностью. 2.7.1. Метод продолжения по ограниченному пара- метру. Выберем параметр rs(rmin, Гтах). Алгоритм (2.58) в этом случае принимает вид W*+i. r^JAv^+^^-I^^, r^,), (2.65) где р^О, 1,...—индекс ньютоновских итераций; k= ==0, 1,...—номер очередного продолжения; r^;+l==/'t+ +Л/'&; (pfe+i—искомая переменная. Значение q'*01 прог- нозируется. В простейшем случае ф^ ==<р*, в более 6E/a, равное номиналь- ному E(t)=E. Алгоритм расчета статического режима совпадает с алгоритмом расчета переходных процессов и имеет вид (2.66) где n=0, I...—номер шага расчета; р—индекс нью- тоновских итераций. Критерием окончания расчета является выполнение условия ||ф<г+1—<рк!1-<е, при котором переходный про- цесс затухает и потенциалы принимают почти постоян- 68 ные значения, соответствующие статическому режиму. Поскольку точность расчета переходного процесса в данном случае не имеет значения, для ускорения расче- та шаг \t==tn+\—tn в (2.66) следует выбирать макси- мально большим, при котором еще сохраняется устой- чивость расчета переходного процесса. С этой целью в (2.66) целесообразно использовать Л-устойчивый не- явный метод Эйлера. Нетрудно видеть, что метод уста- новления (2.66) можно рассматривать как метод про- должения по неограниченному параметру t с использо- ванием движущейся области сходимости. Основные преимущества метода установления перед другими методами заключаются в следующем. Во-первых, алгоритм расчета статического режима (2.66) мало отличается (за исключением выбора шага и критерия окончания расчета) от алгоритма расчета переходного процесса. Благодаря этому сокращается объем и время разработки программы расчета статики и динамики. "Во-вторых, метод установления для большинства схем допускает достаточно произвольный выбор началь- ных условий (ро, поскольку в конечном счете переходный процесс независимо от (ро стремится к одним и тем же асимптотическим значениям. Исключение составляют лишь многоустойчивые схемы, статическое состояние которых зависит от выбора (ро. В-третьих, как правило, расчет статического режи- ма выполняется не только для расчета статических па- раметров схемы, но и весьма часто с целью определения начальных условий для переходных процессов. Метод установления позволяет без дополнительных преобразо- ваний схемы определить начальные значения ее потен- циалов с учетом всех реактивных элементов. Это важно отметить, поскольку, например, статические режимы схемы без индуктивностей и с индуктивностями могут быть различны. Недостатком метода установления следует считать довольно значительное время расчета статического ре- жима. Для его уменьшения используется описываемый ниже метод. 2.7.3. Комбинированный метод. Комбинированный метод продолжения по ограниченному и неограниченно- му параметрам [36] состоит в использовании одновре- менно метода (2.66) и метода (2.65), так что алгоритм можно записать в виде 69 Y(4^i,^, ^^^,=--1(4)^,^, /-^,), (16?) П+1 П+1 fl+1 Реализация алгоритма (2.67) заключается в том, что сначала при г=Гшах производится большой шаг Д< по времени (изменение индекса п). При этом проводимость емкостных ветвей yc=Cf'At становится весьма малой, что соответствует их разрыву, а проводи- мость индуктивных ветвей yL^^t/L—большой, что соответствует режиму, близкому к короткому замыканию. Если ньютоновские итерации (индекс р) при этом сходятся, то расчет окончен, если же нет, то далее алгоритм (2.67) реализуется подобно алгоритму (2.65) по индексу k от /-i==?-min до значений '"t+i='"max, после чего вновь выполняется большой шаг Л< и реали- зуется (2.66). Комбинация (2.65) и (2.66), заключенная в (2.67), повышает надежность программы расчета статического режима, 2.8. Алгоритмы расчета переходных процессов 2.8.1. Алгоритм на основе неявного метода. Расчет переходных процессов в базисе узловых потенциалов выполняется по алгоритму (2.66). При этом использу- ются модели реактивных элементов вида (2.63). Начальные условия (ро определяются из расчета ста- тического режима, алгоритмы автоматического выбора шага, прогноза и контроля локальной точности были описаны в § 2.2. Методика расчета переходных процессов по (2.66) заключается в следующем. Пусть используются про- стейшие дискретизации уравнений (2.63) реактивных ветвей, соответствующие неявному методу Эйлера. По- лагая п==0, р==0, (p^i^cpo и решая систему линейных уравнений Y(q)(())l)ЛcpWl=—Цср^), определяем Atp^i, откуда находим (p(l)l==Лq)(o)l—(p(°)i. Затем вновь решая систему линейных уравнений Y^'^Aqi^i^—HtpOl), определяем Aq)0^, (p^i и т. д. до тех пор, пока выпол- нится какой-либо критерий сходимости ньютоновских итераций. После этого для п=\ выполняется прогноз значений (p^z и снова реализуется метод Ньютона. Когда будет получено достаточное количество «разгонных» точек, можно перейти к расчету переходных процессов мето- дами высокого порядка, например методом ФДН (см. § 2.2) с использование уравнений реактивных элемен- тов, соответствующих дискретизации (2.51). Для прогноза используются (2.43), (2.44) а в случае, метода ФДН — (2.53) при х==ср. 2.8.2. Комбинированный алгоритм. Обычно считает- ся, что использование дискретизации производных пе- 70 ред формированием модели схемы автоматически при- водит к реализации неявных методов численного реше- ния дифференциальных уравнений. Однако если исполь- зовать дискретизацию производных, соответствующую (2.64), то можно реализовать и явные методы [41]. Пусть имеем активную /?С-схему. Разобьем множе- ство всех узлов схемы на подмножества емкостных С, гибридных Г и резистивных R узлов, которым инцидент- ны только емкостные, емкостные и резистивные (в том числе источники) линейные и нелинейные ветви и толь- ко резистивные ветви соответственно. Тогда, исполь- зуя для формирования модели схемы дискретизирован- ные уравнения емкостных ветвей, приведенные к виду (2.64), С "С.П+-1 „кон Уп+i с ,, /, нач „кон^ | г :=Уc(^,-V„+,)+/c,, получаем алгоритм расчета переходных процессов, в котором реализуется какой-либо (в зависимости от спо- соба дискретизации) явный метод решения ОДУ: ffг,Clг+\ L ^".Rn (2.68) Здесь Yc и 1г,с — емкостная матрица узловых проводи- мостей и резистивно-емкостной вектор узловых токов для узлов типа Г и С, образованные всеми емкостными и резистивными ветвями, включенными между узлами типа Г, С или Г и С; Уд, 1д — резистивные матрицы уз- ловых проводимостей и вектор узловых токов для узлов типа R, образованные резистивными ветвями, включен- ными между узлами типа R или R и Г. Таким образом, нижняя подсистема относится толь- ко к чисто резистивным узлам и решается относительно <ряп методом Ньютона, причем в процессе решения по- тенциалы (fm узлов типа Г считаются постоянными. Найденные значения (рдп подставляются в верхнюю под- систему для расчета токов 1г и потенциалов ут+', tpcn+i. Метод (2.68) —явный, так как токи резистивных ветвей вычисляются для момента tn. Следовательно, и токи t'c также будут соответствовать моменту fn, а не tn+\, что и требуется в явных методах. 7} Рис. 2.16. К расчету схем явным методом решения ОДУ В качестве иллюстрации метода (2.68) приведем вид уравнений этого метода для схемы, показанной на рис. 2.16. В этой схеме узлы /, 2, 4—типа Г, узел 3—типа R. Уравнения явного метода узло- вых потенциалов имеют вид Сй<С с,Atс,+с, о о оп Vl.n+i м м с, " о Тг,п+1 == о о м д1 ?<. n+l о о(?!- о-?2)п(с\ дч (С,\ ы, 1 ^ 'P^n+l ^)-/(<,Са г)1 (fa—'hL —(?1—Уг)п [^(«рз +^)——?4)n +?2,П ^С, + Rs К, +У4,п Д< (?s - ч".) •п , ^з.п (^-vJ R. -г ^3 +'(уРз.,>)+ ^ Значения потенциалов Фз,е, ф4,о должны быть заданы как начальные условия. Эти значения подставляются в верхнюю под- систему. Из решения верхней подсистемы определяются (pi,i, <р2,ь (p4,i, которые затем используются как постоянные величины при ре- шении методом Ньютона нижней резистивной подсистемы, состоя- щей в данном случае из одного уравнения для определения q>g,i. Далее процесс решения аналогичен. Легко показать, что данные уравнения соответствуют явному методу Эйлера. Для более сложных явных методов, учитывающих предысторию процесса в точках fn-i, ..., tn-k, в векторе узловых токов появятся дополнительные составляющие в виде источников тока /((pn-i, ..., tpit-ii) для каждой емкостной ветви, 73 Преимущество алгоритма (2.68) перед обычным ме- тодом переменных состояния, также ориентированным на явные методы решения ОДУ, заключается в просто- те формирования уравнений, основанного на классифи- кации узлов при просмотре исходной таблицы, описы- вающей топологию и элементы схемы. Преимущество {2.68) перед алгоритмом Ньютона (2.66) состоит в мень- ших затратах времени на одну итерацию, так как мат- рица Якоби в (2.68) вычисляется только для R-узло'в, количество которых обычно невелико. Если /?-узлы отсутствуют, т. е. к каждому узлу под- ключена хотя бы одна емкость, что характерно, напри- мер, для схем на МДП-транзисторах, то выигрыш (2.68) перед (2.71) в затратах времени и памяти на один шаг расчета переходных процессов становится еще больше, так как отпадает необходимость в матрице Якоби, а верхняя гибридно-емкостная подсистема на каждом временном шаге решается только один раз. Недостатком (2.68) является ограничение на значе- ние шага интегрирования ОДУ, определяемое появле- нием неустойчивости решения и характерное для всех явных методов. Если это ограничение существенно за- медляет расчет, от (2.68) легко перейти к .(2.66), до- формировывая Yc, Уд до полной матрицы схемы Y и пе- реходя к решению методом Ньютона для всей системы уравнений. Следовательно, комбинируя различным об- разом шаги расчета, выполняемые по (2.66) и (2.68), можно добиться существенного сокращения общего вре- мени решения. 2.9. Алгоритм расчета частотных характеристик Метод узловых потенциалов позволяет формировать узловые уравнения не только для временной, но и для частотной области. В этом случае вся предыдущая ме- тодика формирования узловых уравнений сохраняется, изменяются лишь компонентные уравнения реактивных ветвей, принимающие следующий вид: ,J_ h)L (у" (2.69) где ср"34, ср"011—потенциалы на концах реактивных вет- вей; j — мнимая единица; со — частота. 73 Соответственно проводимости реактивных ветвей t/c==JcoC; y^=-j/((oL). (2.70) Уравнения (2.69) используются при формировании вектора узловых токов, а уравнения (2.70)—матрицы узловых проводимостей. При этом постоянные источни- ки напряжения Е закорачиваются, а постоянные источ- ники тока /—размыкаются. В результате получим узло- вое уравнение линейной схемы в частотной области Y(j(o)(p(jco)=-l(io3). (2.71) В отличие от (2.66) для временной области, которое в каждый момент времени нужно решать несколько раз до сходимости, уравнение (2.71) на каждой частоте oi нужно решать лишь один раз, поскольку схема линей- ная. Решение уравнения (2.71) можно выполнить, ис- пользуя стандартную программу либо обращения ком- плексной матрицы, либо решения системы линейных уравнений с комплексными коэффициентами. Подстав- ляя в (2.71) разные значения частоты со, полагая, что комплексная амплитуда входного сигнала равна едини- це, и вычисляя вектор