Введение При выращивании монокристаллов лейкосапфира методом Киропулоса типовой технологический процесс выращивания проходит несколько этапов: получение расплава (разогрев), затравление, рост и охлаждение. Затравление и начало роста являются  самыми ответственными этапами. Здесь часто возникают ошибки, связанные либо с очень большим сбросом мощности оператором, ведущим к залипанию кристалла, либо с очень малым, ведущим к остановке роста и даже отрыву кристалла от расплава [1–3]. Конструкция ростовой установки и метод выращивания не дают возможности непосредственного измерения температуры. Поэтому на первый план выступают косвенные методы, например по изменению показаний измерителя мощности. По характеру изменения мощности на нагревателе можно определить закон изменения температуры [4]. То есть появляется возможность по результатам косвенных измерений управлять процессом выращивания, например используя модель в составе беспоисковой адаптивной системы стабилизации скорости кристаллизации [5]. Для автоматизации данного процесса решаются различные задачи, одной из которых является получение динамических характеристик объекта управления. Причем идентификация в общем случае может быть проведена для неизвестной формы описания модели объекта [6; 7]. В нашем случае с учетом тепловой природы процесса примем, что модель расплава представляет собой апериодическое звено второго порядка. Проведем идентификацию постоянных времени объекта с помощью пакета ident системы инженерных расчетов Matlab [8] по данным напряжения и мощности нагревателя. Применение  способа идентификации Process Models в некоторых случаях даёт неудовлетворительные результаты, поэтому было приято решение использовать для идентификации аппарат пространства состояний, а также проверить влияние начальных условий (НУ) на результаты идентификации.  Математическое описание пространства состояний            В современных пакетах математической обработки форма пространства состояний базируется на матричном представлении  и представляет собой наиболее удобную и компактную форму описания многозонных объектов. В пространстве состояний  математическое описание системы задается в векторно-матричной форме [9]: ,                               (1)где X – вектор состояния размера n×1; U – вектор входных воздействий размера m×1;A – матрица системы размера n×n; B – матрица входа размера n×m.В то же времят  данная форма дает лишь универсальный механизм, конкретное содержание которого зависит от природы описываемого объекта. Рассмотрим представление в форме пространства состояний объекта второго порядка  с передаточной функцией ,          (2)где t1, t2 – постоянные времени объекта;  Kp – коэффициент передачи объекта; Tp1 = t1t2;  Tp2 = t1+t2.Такому объекту соответствует дифференциальное уравнение вида .                     (3)Поскольку кроме текущего состояния объекта необходимо учитывать и скорость его изменения, примем вектор состояния в виде , где P – мощность, выделяемая на нагревателе, Вт; – скорость ее изменения, Вт/с.Тогда  уравнение для пространства состояний  примет вид .                     4)Определим размерности матриц, входящих в это выражение. Вектор-столбец в левой части имеет размерность 2×1, т.е. для сохранения размерности матрица A должна иметь размерность 2×2. Поскольку входное воздействие (напряжение на нагревателе) является одномерной величиной, то матрица B должна иметь размерность 2×1. Выразим эти матрицы в общем виде через их элементы с учетом соглашения системы Matlab [10] (первый индекс кодирует номер строки, а второй – номер столбца): , .С учетом принятых обозначений представим (4) в виде .           (5)Матричное уравнение (5) можно преобразовать к системе дифференциальных уравнений в соответствии с правилами перемножения матриц:                            (6)Определим значения коэффициентов матриц A и  B таким образом, чтобы получить структуру уравнения (3). Поскольку уравнение (3) второго порядка, то в системе (6) нужно оставить только второе уравнение. Для этого первое уравнение в ней сведем к тождеству, т.е. положим a11 = 0, a12 = 1, b1 = 0. Уравнение (2) приведем к виду ,а также перенесем в правую часть все, кроме выражения для второй производной .Тогда возможно определить недостающие коэффициенты матриц A и  B. В итоге они примут вид   , , , D = 0.Идентификация с помощью пространства состояний, влияние НУ Прежде всего проверим адекватность полученной модели. Для этого сравним реакцию модели объекта второго порядка (2), созданную с помощью функции tf, с реакцией модели в форме пространства состояний (State Space). На вход того и другого объекта подадим линейно нарастающий сигнал с помощью функции lsim. Зададимся постоянными времени t1 = 2500 с, t2 = 50 с, коэффициентом передачи Kp = 5. Для случая нулевых НУ графики переходных процессов совпадают (рис. 1).Если же задать ненулевую начальную скорость, например 0,01 Вт/с, то начальная стадия переходного процесса существенно изменится (рис. 1). Также меняется вид графика скорости изменения мощности (рис. 2). В системе Matlab начальные условия задаются с помощью функции initial.После экспорта в рабочее пространство Workspace и преобразования вида tf(n4s2) будет получена передаточная функция .После приведения ее к каноническому виду (1 при S в нулевой степени в знаменателе) получим, чтоKp = 5, Tp1 = 124790, Tp2 = 2548.Учитывая, что t1 = Tp1/ Tp2, t2 = Tp2 - t1, можно утверждать, чтоt1 = 48,978 c, t2 = 2499 c.          Аналогичным образом проведем идентификацию с ненулевыми и нулевыми НУ. Результаты приведены в табл. 1, где FTE – процент соответствия полученной модели исходным данным. Таблица 1Результаты идентификацииНУKpt1, сt2, сFTE, %Нулевые 548,9782499100Ненулевые 548,9782499100  По данным табл. 1 можно сделать вывод, что наличие НУ не влияет на результаты идентификации параметров объекта. Таким образом, при идентификации реального объекта НУ можно не учитывать.  Идентификация экспериментальных данных                      На рис. 4 показан фрагмент экспериментальных данных для напряжения и мощности нагревателя, полученных при выращивании монокристалла. Характерной особенностью данных зависимостей является то, что мощность меняет своё значение с 42,5 до 42,8 кВт, в то время как напряжение поддерживается на уровне 8,15 В (этому предшествовал линейный сброс напряжения). Данный эффект можно объяснить большой инерционностью конвективных потоков в расплаве, которые выносят тепло со дна тигля на поверхность расплава, изменяя распределение температур по высоте вольфрамового высокотемпературного нагревателя. Процесс выращивания монокристалла чувствителен к точности поддержания температуры, и отклонение её значения от поддерживаемой на величину более 0,1 °C является критичным. Поскольку из-за высокой температуры расплава её измерение напрямую невозможно, то используется косвенный метод на основе данных мощности, изменение которой примерно на 60 Вт соответствует 1°C.Кроме приведенного на рис. 4, рассматриваемый процесс выращивания монокристалла  содержит 5 характерных участков со сходной зависимостью мощности от напряжения. Для оценки стабильности методов идентификации Process Models и Space State проведем обработку всех пяти участков, а полученные данные сведём в табл. 2 и 3 соответственно. При анализе результатов будем учитывать не только  процент соответствия полученной модели исходным данным (FTE), но и разброс получаемых значений, а также ожидаемый исходя из теплофизических характеристик ростовой установки диапазон постоянных времени (для t1 – десятки секунд, для t2 – тысячи секунд). Таблица 2Результаты идентификации с помощью Process ModelsУчастокKp FTE, %12,12277,1902290696,6623,220,007701198,4835,220,91810-636,1145,110,028164074,2854,934,697116941,42    Таблица 3Результаты идентификации с помощью State SpaceУчастокKp FTE, %12,146,7722344598,0223,2337,408678197,7535,2352,574149054,4945,1277,175136271,6054,9378,647143840,23 Заключение Сравнение значений t2 в табл. 2 и 3 показывает, что идентификация при помощи Space State позволяет получить более достоверные результаты, о чём свидетельствует отсутствие значений, резко отличающихся от значений на соседних участках. Также получены более корректные значения t1, поскольку предполагается, что она равна десяткам секунд. Полученная модель позволяет косвенно оценить интенсивность движения тепловых потоков, создающих конвективный теплообмен, что может быть использовано в качестве корректирующего воздействия в системе управления установкой выращивания монокристаллов.