Two-channel time-optimal control of induction heating process with maximum temperature constraint

Full Text

Процесс нагрева металлических изделий и полуфабрикатов уже давно является неотъемлемой частью производственного процесса в различных отраслях промышленности. Технологии индукционного нагрева используются для подготовки металлических заготовок к пластической деформации, для чего их сначала предварительно нагревают до требуемой температуры [1].

Нередки случаи, когда возникает необходимость в коррекции оптимальных режимов нагрева, вызванной жесткими технологическими требованиями к поведению температурного поля в процессе управления. Подобные ограничения могут накладываться на максимальное значение температуры по всему объему нагреваемой заготовки во избежание нежелательных изменений структурных свойств металла.

В рамках данной работы на первом этапе будет рассмотрена такая постановка задачи оптимального управления (ЗОУ), при которой фазовое ограничение не будет учитываться, чтобы определить, выполняется ли в таком варианте требование к максимальной температуре любой точки заготовки «автоматически». Если оно не выполняется, то необходимо проводить коррекцию оптимального режима индукционного нагрева с учетом фазового ограничения на температурное поле.

Как известно, решение задач, связанных с процедурой предварительного нагрева заготовки до требуемого температурного состояния, в большинстве своем требует привлечения теории оптимального управления объектами с распределенными параметрами (ОРП).

Общая постановка задачи оптимального управления

В качестве объекта управления рассматривается технологический процесс индукционного нагрева металлической заготовки с двумя управлениями по мощности внутренних источников тепла, воздействующих на различные поверхности неограниченной пластины. В условиях пренебрежения неравномерностью распределения температуры по длине и ширине заготовки температурное поле пластины  описывается в зависимости от пространственной координаты x и времени t линейным пространственно-одномерным однородным уравнением теплопроводности следующего вида [2]:

$\frac{\partial Q(x,t)}{\partial t}=a\frac{{\partial }^{2}Q(x,t)}{\partial x^2}+W_1\left(x\right)u_1\left(t\right)+W_2\left(x\right)u_2\left(t\right),x\in \left[\mathrm{0,}R\right];t\in \left[\mathrm{0,}T\right].$ (1)

Здесь a – коэффициент температуропроводности; R – толщина пластины; $W_1\left(x\right),W_2\left(x\right)$ – функция распределения электромагнитных внутренних источников тепла; $u_1\left(t\right),u_2\left(t\right)$ – удельные мощности тепловыделения, рассматриваемые в качестве управляющих воздействий и изменяющиеся в заранее заданных пределах:

$0\le u_1\left(t\right)\le u_{1\max };0\le u_2\left(t\right)\le u_{2\max }.$ (2)

Начальное распределение принимается равномерным по всему объему пластины и, в частности, равным нулю:

$Q\left(x\mathrm{,0}\right)=Q_0\left(x\right)=Q_0=\mathrm{const}=0.$(3)

На обеих поверхностях действуют одинаковые мощности внутренних источников тепла $W_1\left(x\right)$ и $W_2\left(x\right)$, определяемые известными выражениями [3]:

$W_1\left(x\right)=\frac{\text{ch}\left(\sqrt{2}\xi \frac{x}{R}\right)-\cos \left(\sqrt{2}\xi \frac{x}{R}\right)}{\text{sh}\left(\sqrt{2}\xi \right)-\text{sin}\left(\sqrt{2}\xi \right)}\sqrt{2}\xi ,W_2\left(x\right)=\frac{\text{ch}\left(\sqrt{2}\xi \frac{R-x}{R}\right)-\cos \left(\sqrt{2}\xi \frac{R-x}{R}\right)}{\text{sh}\left(\sqrt{2}\xi \right)-\text{sin}\left(\sqrt{2}\xi \right)}\sqrt{2}\xi ,$

где ξ – характерный параметр, определяемый соотношением

$\xi =\frac{R\sqrt{2}}{\delta },\delta =\sqrt{\frac{2}{\mu \omega \sigma }}.$

Здесь δ – глубина проникновения тока в металл; ω – частота питающего тока; σ – электропроводность нагреваемого материала; μ – абсолютная магнитная проницаемость. Типовые граничные условия для модели объекта (1) - (2) в общем случае имеют вид [4]

$\begin{array}{c}-\lambda \frac{\partial Q\left(\mathrm{0,}t\right)}{\partial x}=Q_{ср}\left(t\right)-{\alpha }_{1}Q\left(\mathrm{0,}t\right),t>\text{0;}\\ \lambda \frac{\partial Q(R,t)}{\partial x}=Q_{ср}\left(t\right)-{\alpha }_{2}Q(R,t),t>\text{0,}\end{array}$(4)

где    $\lambda$   – коэффициент теплопроводности;

${\alpha }_1,{\alpha }_2$ – заданные теплофизические постоянные;

$Q_{ср}\left(t\right)=0$ – температура окружающей среды.

Соотношения (4) соответствуют граничным условиям 3-го рода [2]. Наглядное представление управляющих воздействий на поверхности заготовки представлено на рис. 1.

 

Рис. 1. Иллюстрация двухканального распределенного управления

 

В момент T окончания процесса управления требуется обеспечить заданную точность ε равномерного приближения конечного распределения температуры $Q(x,t)$ к заданному $Q^{*}\left(x\right)=Q*=const>Q_0$:

$\underset{x\in \left[\mathrm{0,}R\right]}{\max }\left|Q(x,T)-Q*\right|\le \epsilon .$ (5)

На всем протяжении периодического нагрева может накладываться ограничение на максимальное значение температуры по всему объему заготовки, нарушение которого вызывает оплавление нагреваемой заготовки или приводит к недопустимому изменению структуры металла. Такое ограничение заключается в требовании к максимальной температуре не превысить некоторую заданную допустимую величину $Q_{\text{доп}}$ и описывается в форме следующего неравенства:

$Q\left(x\right)=\max Q(x,t)\le Q,t\in \left[\mathrm{0,}T\right].$ (6)

Как уже было сказано, на первом этапе решения ЗОУ требование (6) не будет приниматься в расчет при постановке задачи.

В качестве критерия оптимальности выступает общее время T процесса нагрева в виде следующего интегрального функционала качества:

$I=\underset{0}{\overset{T}{\int }}dt=T\to \underset{u_1\left(t\right),u_2\left(t\right)}{\min }.$ (7)

Объект, описываемый соотношениями (1) - (7), является объектом с распределенными параметрами. Для поиска вход-выходной характеристики ОРП с двумя управляющими воздействиями применяется метод конечных интегральных преобразований. Применение этого метода, позволяющего описать ОРП соответствующим модальным представлением, подробно описано в [5] при решении подобной задачи.

Метод конечных интегральных преобразований [4-7] определяет функцию состояния температурного поля $Q(x,t)$ в зависимости от пространственной координаты $x\in \left[\mathrm{0,}R\right]$ и времени $t\in \left[\mathrm{0,}T\right]$ в форме разложения в ряд по собственным функциям ${\overline{Q}}_n({\mu }_n,t)$ с временными модами ${\overline{Q}}_n({\mu }_n,t)$ начально-краевой задачи (1) - (5):

$\frac{d{\overline{Q}}_n\left({\mu }_n,t\right)}{dt}=-{{\mu }_n}^2{\overline{Q}}_n\left({\mu }_n,t\right)+\frac{1}{\lambda E_n}\left(W_{1n}{u}_1\left(t\right)+W_{2n}{u}_2\left(t\right)\right)\text{, }n=\mathrm{1,2,}\mathrm{...};$ (8)

${\overline{Q}}_n\left({\mu }_n\mathrm{,0}\right)={\overline{Q}}_0\left({\mu }_n\right)=0;$ (9)

$Q(x,t)=\sum _{n=0}^{\infty }{\overline{Q}}_n({\mu }_n,t)\phi ({\mu }_n,x),$ (10)

где собственные функции $\phi ({\mu }_n,x)$, собственные числа ${\mu }_n^2$, нормирующие множители $\frac{1}{E_n}$ и числа ${\eta }_n$ определяются следующими соотношениями [4, 5]:

$\begin{array}{l}{\phi }_n({\mu }_n,x)=\frac{1}{E_n}\cos \left({\eta }_n\frac{x}{R}\right),n=\mathrm{1,2,...};\\ {\eta }_n\text{\hspace{0.17em}tg}{\eta }_n-Bi=0;{\mu }_n=\frac{\sqrt{a}}{R}{\eta }_n;\\ E_n=\sqrt{\frac{a}{R}\cdot \frac{2{\eta }_n}{{\eta }_n+\sin {\eta }_n\cos {\eta }_n}},n=\mathrm{1,2,...}\end{array}$

Выражение для $Q(x,t)$ в форме (10) находится в явном виде по алгоритму, описанному в [5].

Зависимость температурного поля от управляющих воздействий по мощности внутренних источников тепла в выражении (10) без учета фазовых ограничений описывается в виде следующего равенства:

$\begin{array}{c}Q(x,t)=\sum _{n=1}^{\infty }{\overline{Q}}_n({\eta }_n,t)\phi ({\eta }_n,x)=\sum _{n=1}^{\infty }\phi ({\eta }_n,x)\times \\ \times \left.\left[\underset{0}{\overset{R}{\int }}\underset{0}{\overset{t}{\int }}\left(W_1\left(\xi \right)\phi ({\eta }_n,\xi )\frac{1}{a}{u}_1\left(\tau \right)+W_2\left(\xi \right)\phi ({\eta }_n,\xi )\frac{1}{a}{u}_2\left(\tau \right)\right)\text{exp}\left(-{\eta }_n^2\frac{a}{R^2}(t-\tau )\right)d\xi d\tau \right.\right].\end{array}$ (11)

Выражение (11) представлено в форме двойных интегралов по пространственно-временной координате и определяет реакцию объекта на воздействия $W_1\left(x\right)u_1\left(t\right),W_2\left(x\right)u_2\left(t\right).$

В итоге может быть сформулирована следующая задача оптимального по быстродействию управления. Требуется найти такие стесненные ограничением (2) программные управляющие воздействия $u_1^{*}\left(t\right),u_2^{*}\left(t\right)$, которые переводят объект управления (8) из заданного начального состояния (9) в требуемое конечное согласно (5), где $Q(x,t)$ определяется выражением (11), при минимальном значении критерия оптимальности (7) в условиях, когда ограничение на поведение температурного поля (6) не учитывается на всем протяжении процесса управления.

Параметризация управляющих воздействий

На основании процедуры принципа максимума Понтрягина [8], применение которого распространяется на бесконечномерный объект управления [9] (9), (10), оптимальное по быстродействию управление без учета фазовых ограничений устанавливается в форме релейных функций времени, попеременно принимающих только свои предельно допустимые значения согласно (2) с точностью до числа $N_1$ и $N_2$ и длительностей ${\Delta }_{1i}^{\left(N_1\right)}$ и ${\Delta }_{2i}^{\left(N_2\right)}$ интервалов постоянства управляющих воздействий [3-9] соответственно $u_1\left(t\right)$ и $u_2\left(t\right)$:

$\begin{array}{c}{u}_1^{*}({\Delta }_1^{\left(N_1\right)},t=\frac{u_{1\max }}{2}+{\left(-1\right)}^{j+1}\frac{u_{1\max }}{2},\sum _{i=0}^{\text{j-1}}{\Delta }_{1i}^{\left(N_1\right)}<t<\sum _{i=0}^{\text{j}}{\Delta }_{1i}^{\left(N_1\right)},j=\overline{\mathrm{1,}{N}_1}\text{, }{\Delta }_{10}=\mathrm{0,}{\Delta }_1^{\left(N_1\right)}={\Delta }_{1i}^{\left(N_1\right)};\\ u_2^{*}({\Delta }_2^{\left(N_2\right)},t)=\frac{u_{2\max }}{2}+{\left(-1\right)}^{j+1}\frac{u_{2\max }}{2},\sum _{i=0}^{\text{j-1}}{\Delta }_{2i}^{\left(N_2\right)}<t<\sum _{i=0}^{\text{j}}{\Delta }_{2i}^{\left(N_2\right)},j=\overline{\mathrm{1,}{N}_2},{\Delta }_{20}=\mathrm{0,}{\Delta }_2^{\left(N_2\right)}={\Delta }_{2i}^{\left(N_2\right)}.\end{array}$ (12)

В пределах первого интервала постоянства согласно ограничениям (2) здесь и далее принимается $u_1^{*}\left(t\right)=u_{1\max },u_2^{*}\left(t\right)=u_{2\max }$ в рассматриваемой задаче нагрева пластины до температуры $Q^{*}>Q_0=0.$

Выражение (12) свидетельствует о том, что в процессе нагрева объекта происходит скачкообразное (релейное) переключение управляющих воздействий с максимального значения, соответствующего нагреву заготовки с максимальной интенсивностью, когда каждый из индукторов включается на полную мощность, на минимальное значение мощности при отключении индукторов от сети.

Условимся, что далее будет рассмотрен наиболее употребительный на практике случай двухинтервального управления (рис. 2б), для этого примем $N_1=N_2=2$. Поскольку оба управляющих воздействия (12) должны заканчиваться в один и тот же момент времени, то длительности интервалов каждого из управлений стесняются условием равенства отдельных сумм: ${\Delta }_{11}^{\left(2\right)}+{\Delta }_{12}^{\left(2\right)}={\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)}$. После подстановки управления вида (12) в (11) и вычисления интегралов при $t=T$ получим параметризованное представление конечного температурного состояния:

$\begin{array}{c}Q(x,{\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)})=\sum _{n=1}^{\infty }\frac{\lambda R^2}{a{E}_n{\eta }_n^2}\left(\cos \left({\eta }_n\frac{x}{R}\right)+\frac{B{i}_1}{{\eta }_n}\sin \left({\eta }_n\frac{x}{R}\right)\right)\times \\ \times \left[W_{1n}{u}_{1\max }\left(e^{-{\eta }_n^2\frac{a}{R^2}{\Delta }_{12}^{\left(2\right)}}-e^{-{\eta }_n^2\frac{a}{R^2}({\Delta }_{11}^{\left(2\right)}+{\Delta }_{12}^{\left(2\right)})}\right)\right.+\left.W_{2n}{u}_{2\max }\left(e^{-{\eta }_n^2\frac{a}{R^2}{\Delta }_{22}^{\left(2\right)}}-e^{-{\eta }_n^2\frac{a}{R^2}({\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)})}\right)\right],\end{array}$ (13)

где $W_{1n},W_{2n}$ – моды разложения функций $W_1\left(x\right),W_2\left(x\right)$ в ряд по собственным функциям.

Релейный характер двухканального двухинтервального оптимального по быстродействию управления, найденный в (12), представлен на рис. 2а, из которого видно, что действие обоих управляющих воздействий начинается и заканчивается одномоментно. В связи с этим удобнее перейти к двум комбинациям искомых величин: в качестве трех независимо варьируемых параметров выберем ${\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}$ и ${\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)},$ а интервалы нагрева минимальной интенсивности для первого ${\Delta }_{12}^{\left(2\right)}$ и второго ${\Delta }_{22}^{\left(2\right)}$ управления обозначим в (13) как ${\Delta }_{22}^{\left(2\right)}={\Delta }_{11}^{\left(2\right)}+{\Delta }_{12}^{\left(2\right)}-{\Delta }_{21}^{\left(2\right)}$ и ${\Delta }_{12}^{\left(2\right)}={\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)}-{\Delta }_{11}^{\left(2\right)}$.

 

Рис. 2. Характер изменения оптимальных управляющих воздействий (а) и управляемой величины (б) в процессе индукционного нагрева

 

Редукция к задаче полубесконечной оптимизации

При полученном параметрическом представлении искомых управляющих воздействий (13) критерий оптимальности (7) в соответствии с рис. 2а определяется в виде простой суммы длительностей отдельных интервалов постоянства оптимального управления:

$I={\Delta }_{11}^{\left(2\right)}+{\Delta }_{12}^{\left(2\right)}={\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)}\to \min ,$ (14)

а условие (5) оценки конечного распределения температур будет иметь вид:

$\Phi \left({\Delta }_{ij}^{\left(2\right)}\right)=\underset{x\in \left[\mathrm{0,}R\right]}{\max }\left|Q\left({\Delta }_{ij}^{\left(2\right)}\right)-Q*\right|\le \epsilon ,i,j=\overline{\mathrm{1,2}},\epsilon >0.$ (15)

Здесь $Q(x,{\Delta }_{ij}^{\left(2\right)})$ определяется по формуле (13).

Таким образом, производится точная редукция исходной ЗОУ к задаче полубесконечной оптимизации (ЗПО) на минимум целевой функции (14) конечного числа переменных ${\Delta }_{ij}^{\left(2\right)},i,j=\overline{\mathrm{1,2}}$ с бесконечным числом ограничений (15), порождаемых требованием обеспечения заданной точности равномерного приближения во всех точках $x\in \left[\mathrm{0,}R\right]$ [9, 10].

Решение задачи полубесконечной оптимизации по схеме альтернансного метода

Задача (14), (15) разрешима в условиях $N_1=N_2=2$ не при всех , а только для $\epsilon \ge {\epsilon }_{\min }^{\left(2\right)}$ в (15), где ${\epsilon }_{\min }^{\left(2\right)}$ – минимально достижимая величина  в рассматриваемом классе граничных управлений:

${\text{ε}}_{\text{min}}^{\text{(2)}}=\underset{{\Delta }_{ij}^{\left(2\right)}}{\min }\left\{\underset{x\in \left[\mathrm{0,}R\right]}{\max }\left|Q\left({\Delta }_{ij}^{\left(2\right)}\right)-Q*\right|\right\}.$ (16)

Решение ${\Delta }_{ij}^{\left(2\right)},i,j=\overline{\mathrm{1,2}}$ задачи (14), (15) при достаточно малостеснительных допущениях обладает базовыми альтернансными свойствами [11].

Основное свойство заключается в том, что число точек $x_j^0\in \left[\mathrm{0,}R\right],j=\overline{\mathrm{1,}{R}_x}$, в которых достигаются предельно допустимые абсолютные отклонения $Q(x,{\Delta }_{ij}^{\left(2\right)})$ от $Q^{*}$, равные ε, всегда оказывается равным числу всех неизвестных в ЗПО (14), (15).

В качестве неизвестных фигурируют длительности интервалов постоянства ${\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}$ или ${\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}$ оптимального управления и величина минимакса ${\epsilon }_{\min }^{\left(2\right)}$ в случае $\epsilon ={\epsilon }_{\min }^{\left(2\right)}$ в (16).

Таким образом, на искомом решении  задачи быстродействия выполняются  равенств вида

$\left|Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)-Q*\right|=\epsilon ,j=\overline{\mathrm{1,}{R}_x};$ (17)

$R_x=\left\{\begin{array}{l}s,\epsilon >{\epsilon }_{\min }^{\left(2\right)};\\ s+\mathrm{1,}\epsilon ={\epsilon }_{\min }^{\left(2\right)}.\end{array}\right.$ (18)

Здесь s – число свободно варьируемых параметров в составе ${\Delta }_{ij}^{\left(2\right)},i,j=\overline{\mathrm{1,2}},$, равное $s=N_1+N_2-1$ в условиях одинаковой длительности процесса управления для обоих управляющих воздействий (см. рис. 2). В рассматриваемом случае двухинтервального характера управления при $N_1=N_2=2$ имеем $s=N_1+N_2-1=2+2-1=3$ в предельном случае $\epsilon ={\epsilon }_{\min }^{\left(2\right)}$. То есть искомыми параметрами будут $s=3$ длительности интервалов нагрева.

Если требуемая точность температурного отклонения ε задается конкретным числом, превышающим минимаксное значение ${\epsilon }_{\min }^2$, то задача характеризуется только искомыми параметрами в составе ${\Delta }_{ij}^{\left(2\right)},i,j=\overline{\mathrm{1,2}}$, и тогда число точек $R_x$ оказывается равным s согласно условию (18). В этом случае в качестве неизвестных будут фигурировать длительности интервалов постоянства оптимального управления.

Если же значение ε по заданным изначально условиям должно совпадать с заранее неизвестной величиной минимакса ${\epsilon }_{\min }^{\left(2\right)}$, то в соответствии с (18) точек $x_j$ должно быть на одну больше, то есть $R_x=s+1$. Тогда наряду с длительностями интервалов ${\Delta }_{ij}^{\left(2\right)},i,j=\overline{\mathrm{1,2}}$ величину ${\epsilon }_{\min }^{\left(2\right)}$ надо отнести к искомым параметрам.

Рассмотрим наиболее интересный предельный вариант при $\epsilon ={\epsilon }_{\min }^{\left(2\right)}$, тогда число уравнений в (18) будет равно

$R_x=s+1=3+1=\mathrm{4,}$ (19)

при наличии четырех неизвестных: ${\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)},{\epsilon }_{\min }^{\left(2\right)}$ или ${\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)},{\epsilon }_{\min }^{\left(2\right)}.$

Таким образом, соотношения (17) с учетом (18), (19) оказываются замкнутыми относительно всех параметров процесса управления, что указывает на возможность решения исходной задачи оптимального управления.

Основное затруднение теперь состоит в том, что равенствам (17) формально соответствует множество вариантов по форме кривой пространственного распределения $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)$. Для однозначного определения вида этой кривой нужно установить знаки разностей $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)-Q*$ в каждом из уравнений и найти координаты точек $x_j^0$. Эта задача может быть решена только при известной конфигурации кривой температурного распределения $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)-Q*$ на отрезке $\left[\mathrm{0,}R\right]∍x$ при двухинтервальном оптимальном управлении по мощности источников тепла, устанавливаемой на основании физических закономерностей процессов нестационарной теплопроводности в зависимости от величины ε.

Анализ этих закономерностей [11, 12] приводит к двум вариантам конфигурации кривых температурного распределения $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)-Q*$ в условиях $\epsilon ={\epsilon }_{\min }^{\left(2\right)}$ в (17), представленным на рис. 3, где для кривой а принимаются длительности ${\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}$, а для кривой б – длительности ${\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}$.

Отметим, что на кривой а и кривой б, указанных на рис. 3, максимально допустимые отклонения $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)$ от $Q^{*}$ оказываются знакочередующимися в последовательно расположенных точках $x_j\in \left[\mathrm{0,}R\right]$, образуя чебышевский альтернанс.

Теперь, когда характер зависимости $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right)-Q*$ найден из знаний предметной области процесса нестационарной теплопроводности, система (17) может быть однозначно определена для каждого из вариантов кривых конечного температурного распределения.

 

Рис. 3. Два варианта кривых результирующего температурного распределения при $\mathit{\epsilon }\mathbf{=}{\mathit{\epsilon }}_{\mathbf{m}\mathbf{i}\mathbf{n}}^{\mathbf{\left(}\mathbf{2}\mathbf{\right)}}$:

а – при поиске двух интервалов для управления u₁;

б – при поиске двух интервалов для управления u₂

 

Условимся, что будем рассматривать оба варианта кривой температурного распределения, приведенной на рис. 3, для каждой из которых запишем соответствующую систему уравнений, дополняемую условием существования экстремума во внутренних точках отрезка [0, R] максимума или минимума разности температур:

$\begin{array}{l}Q\left(\mathrm{0,}{\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}\right)-Q^{\ast }=-{\epsilon }_{\min }^{\left(2\right)};\\ Q\left(x_2,{\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}\right)-Q^{\ast }={\epsilon }_{\min }^{\left(2\right)};\\ Q\left(x_3,{\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}\right)-Q^{\ast }=-{\epsilon }_{\min }^{\left(2\right)};\\ Q\left(x_4,{\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}\right)-Q^{\ast }={\epsilon }_{\min }^{\left(2\right)};\\ \frac{\partial Q\left(x_j,{\Delta }_{11}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)},{\Delta }_{21}^{\left(2\right)}\right)}{\partial x}=\mathrm{0,}j=\mathrm{2,3,4}\end{array}$ (20, а)

 

$\begin{array}{l}Q\left(x_2,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}\right)-Q^{\ast }={\epsilon }_{\min }^{\left(2\right)};\\ Q\left(x_3,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}\right)-Q^{\ast }=-{\epsilon }_{\min }^{\left(2\right)};\\ Q\left(x_4,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}\right)-Q^{\ast }={\epsilon }_{\min }^{\left(2\right)};\\ Q\left(R,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}\right)-Q^{\ast }=-{\epsilon }_{\min }^{\left(2\right)};\\ \frac{\partial Q\left(x_j,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}\right)}{\partial x}=\mathrm{0,}j=\mathrm{2,3,4}\end{array}$(20, б)

Каждая из систем (20) состоит из семи уравнений с семью неизвестными и решается стандартными численными методами. Численные результаты решения систем уравнений (20, а) и (20, б) с использованием программной среды MATLAB [14-17] были получены применительно к исходным данным, указанным в табл. 1 и 2.

 

Таблица 1  Характеристики нагреваемых заготовок

Параметр	Значение
Материал заготовки	титан
R, толщина заготовки, м	0,2
λ, коэффициент теплопроводности, Вт/(м ∙ ºС)	14
α1, коэффициент конвективной теплопередачи	32,5
α2, коэффициент конвективной теплопередачи	22,5
a, коэффициент температуропроводности, м2/c	4,3∙10-6

 

Таблица 2  Исходные данные для процесса нагрева

Параметр	Значение
Bi1, критерий Био	0,4643
Bi2, критерий Био	0,3214
Q0 начальная температура, ºС	0
Q* конечная температура, ºС	960
P1max максимальная величина поверхностной мощности нагрева по первому каналу управления, кВт/м2	106
P2max максимальная величина поверхностной мощности нагрева по второму каналу управления, кВт/м2	180

 

Результаты расчета оптимального по быстродействию индукционного процесса управления без учета технологических ограничений на максимальную температуру для системы (20, а), соответствующей кривой а на рис. 3: ${\Delta }_{11}^{\left(2\right)}=1524с, {\Delta }_{12}^{\left(2\right)}=606c, {\Delta }_{21}^{\left(2\right)}=1853с, {\epsilon }_{\min }^{\left(2\right)}=\mathrm{16,54}°C.$

Результаты решения системы (20, б), соответствующей кривой б на рис. 3: на рис. 3: ${\Delta }_{11}^{\left(2\right)}=1732с, {\Delta }_{12}^{\left(2\right)}=312с, {\Delta }_{21}^{\left(2\right)}=1416с, {\epsilon }_{\min }^{\left(2\right)}=\mathrm{13,67}°C.$

Полученные при расчетах кривые температурного отклонения от заданной величины представлены на рис.4.

 

Рис. 4. Кривые температурного распределения при $\mathit{\epsilon }\mathbf{=}{\mathit{\epsilon }}_{\mathbf{m}\mathbf{i}\mathbf{n}}^{\mathbf{\left(}\mathbf{2}\mathbf{\right)}}$, полученные при решении систем (20, а) и (20, б) без учета фазового ограничения: а – при поиске двух интервалов для управления u₁; б – при поиске двух интервалов для управления u₂

 

Решение систем уравнений производилось с учетом первых 30 членов бесконечного ряда в выражении (13).

Таким образом, решение системы (20, б) является оптимальным и по быстродействию, и по точности нагрева. Именно для этого варианта решения построим график изменения температуры в процессе нагрева. На рис. 5 можно увидеть характер изменения температурного поля в трех точках: на крайних границах $x=\mathrm{0,}x=R$ и в центре пластины.

Остановимся на рис. 5 подробнее. Можно видеть, что максимальная температура достигает значения более 1500 ºС. Очевидно, что если такой оптимальный процесс с управлением вида (12), рассчитанный без учета (6), не нарушает этого ограничения, то он является оптимальным и с учетом данных ограничений, которые выполняется в некотором роде «автоматически».

Результаты решения задачи оптимального процесса управления во времени без учета технологических ограничений, отвечающие алгоритму (12), указывают, в пределах какого интервала постоянства оптимального управляющего воздействия может возникнуть нарушение ограничения на максимум допустимой температуры и как его следует изменить, чтобы выполнить неравенство (6). То есть первый этап решения задачи заключается в обнаружении выполнения или невыполнения условия (6) в каждом отдельном случае.

 

Рис. 5. Изменение управляемой величины и управляющих воздействий на всем протяжении оптимального процесса управления без учета фазового ограничения

 

Рассмотрим случай, когда условие (6) не выполняется «автоматически» и возникает необходимость ограничения, например, управления $u_2^{*}\left(t\right)$ из-за превышения допустимого уровня температуры $Q_{\text{доп}}=1400$ ºС. Очевидно, что коррекция оптимального режима нагрева в связи с фазовым ограничением оказывается необходимой только в пределах первого интервала процесса нагрева с максимальной мощностью $u_2^{*}\left(t\right)=u_{2\max }$, что соответствует наиболее распространенному на практике случаю.

В источнике [3] показано, что процедура коррекции сводится к подбору таких управляющих воздействий на соответствующих участках в пределах первого интервала, которые поддерживают максимальную температуру $Q_{\max }$ на предельно допустимом уровне $Q_{\text{доп}}$. В [3] предлагается несколько вариантов оптимальной программы мощности нагрева, выбор конкретной из них зависит от специфики отдельной рассматриваемой задачи.

В итоге алгоритм оптимального управления усложняется на первом интервале наличием «участка движения по ограничениям», когда $Q_{\max }=Q_{\text{доп}}$. На втором интервале нагрева в данной постановке задачи оптимальный режим нагрева остается неизменным.

Физика построения нового алгоритма управления с учетом фазового ограничения заключается в следующем. На первом этапе проводится форсированный нагрев заготовки с максимальной интенсивностью при $u_2^{*}\left(t\right)=u_{2\max }$ применительно к рассмотренной выше задаче. В момент, когда температура $Q(x,t)$ достигает величины $Q_{\text{доп}}$, следует участок поддержания на этом уровне управляющим воздействием ${\overline{u}}_2\left(t\right)$ по мощности нагрева, которое скачкообразно уменьшается относительно $u_{2\max }$ с целью снижения скорости роста температуры. Поддержание равенства $Q_{\max }=Q_{\text{доп}}$ требует непрерывного уменьшения мощности тепловыделения по мере выравнивания температурных перепадов.

В соответствии со сказанным алгоритм оптимального по быстродействию управления с учетом фазовых ограничений (6) вместо (12) принимает следующий вид:

$u_2^{*}=\left\{\begin{array}{l}{u}_{2\max },t\in \left(\mathrm{0,}{\phi }_0\right);\\ \overline{u}\left(t\right)\text{, }t\in ({\phi }_0,{\Delta }_{21}^{\left(2\right)});\\ \mathrm{0,}t\in ({\Delta }_{21}^{\left(2\right)},T);\end{array}\right.$ (21)

где $\overline{u}\left(t\right)$ – стабилизирующее управление, которое находится [14]:

$\overline{u}\left(t\right)\cong a_Q+b_Q{e}^{-{\beta }_Q(t-{\upsilon }_0)}.$ (22)

Здесь ${\phi }_0$ – момент достижения $Q_{\text{доп}}$. Коэффициенты $a_Q,b_Q,{\beta }_Q$ находятся по выражениям [14]:

$\begin{array}{c}{a}_Q=\frac{L_1{\mu }_2^2+L_2{\mu }_1^2}{E_1{\mu }_2^2+E_2{\mu }_1^2};b_Q=\frac{\left({\mu }_2^2-{\mu }_1^2\right)\left(E_1{L}_2+E_2{L}_1\right)}{\left(E_1+E_2\right)\left(E_1{\mu }_2^2+E_2{\mu }_1^2\right)};{\beta }_Q=\frac{E_1{\mu }_2^2+E_2{\mu }_1^2}{E_1+E_2};\\ L_n={\theta }_n{A}_n{\eta }_n^2-\frac{u_{1\max }}{u_{2\max }};{\theta }_n=\underset{o}{\overset{R}{\int }}Q\left(\frac{x}{R}{\upsilon }_0\right)\cos \left(\frac{x}{R}{\eta }_n\right)dx;\\ A_n=\frac{2{\eta }_n^2\cos \left({\eta }_{n}R\right)}{\left({\eta }_2^2+B{i}_2^2+B{i}_2\right)\sin \left({\eta }_n\right)};E_n=A_n\cdot W_n\left(\xi \right),\end{array}$ (23)

где ${\eta }_1^2$, ${\eta }_2^2$ – первые два из собственных числе в (10). Именно в выражении (23) для коэффициента $L_n$ заключена особенность корректировки оптимального режима нагрева при двухканальном управлении по сравнению с выражением, представленным в [14], за счет управляющего воздействия $u_{1\max }$ по первому каналу, в котором управление не ограничивается.

Первый временной момент достижения допустимой величины ${\upsilon }_0$ может быть определен по рис. 5 и равен  ${\upsilon }_0=1149$с.

При известных ${\upsilon }_0,\overline{u}\left(t\right)$ результирующее температурное поле вместо (13) вычисляется по следующей формуле для конечного момента времени процесса после подстановки управления вида (21):

$\begin{array}{c}Q(x,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{12}^{\left(2\right)})=\sum _{n=1}^{\infty }\frac{\lambda R^2}{a{E}_n{\eta }_n^2}\left(\cos \left({\eta }_n\frac{x}{R}\right)+\frac{B{i}_1}{{\eta }_n}\sin \left({\eta }_n\frac{x}{R}\right)\right)\times \\ \times \left[W_{1n}{u}_{1\max }\left(e^{-{\eta }_n^2\frac{a}{R^2}{\Delta }_{11}^{\left(2\right)}}-e^{-{\eta }_n^2\frac{a}{R^2}({\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)}-{\Delta }_{11}^{\left(2\right)})}\right)\right.+W_{2n}\left(-1+e^{{\eta }_n^2{\upsilon }_0}+\right.\\ \left.\left.+a_Q\left(e^{{\eta }_n^2{\Delta }_{22}^{\left(2\right)}}-e^{{\eta }_n^2{\upsilon }_0}\right)+\frac{b_Q{\eta }_n^2}{{\eta }_n^2-{\beta }_Q}\left(e^{{\eta }_n^2-{\beta }_Q\left({\Delta }_{22}^{\left(2\right)}-{\upsilon }_0\right)}-e^{{\eta }_n^2{\upsilon }_0}\right)\right)e^{-{\eta }_n^2\left({\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)}\right)}\right].\end{array}$ (24)

В (24) отражен тот факт, что коррекция режима оптимального управления с помощью стабилизирующего управления $\overline{u}\left(t\right)$ накладывается согласно динамике изменения температурного поля (см. рис. 5) на максимальную мощность второго источника тепла $u_{2\max }$ на интервале интенсивного нагрева.

Теперь по описанной ранее схеме альтернансного метода остается решить систему (20, б), где $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right),i,j=\overline{\mathrm{1,2}}$ определяется выражением (24).

 

Рис. 6. Кривая температурного отклонения при $\mathit{\epsilon }\mathbf{=}{\mathit{\epsilon }}_{\mathbf{m}\mathbf{i}\mathbf{n}}^{\mathbf{\left(}\mathbf{2}\mathbf{\right)}}$ полученная при расчетах с учетом фазового ограничения при коррекции управления ${\mathit{u}}_{\mathbf{2}\mathbf{m}\mathbf{a}\mathbf{x}}$

 

Результаты расчета оптимального по быстродействию процесса индукционного нагрева заготовки с ограничением на максимальную температуру при коррекции управления $u_{2\max }$: ${\Delta }_{21}^{\left(2\right)}=1560с, {\Delta }_{22}^{\left(2\right)}=336с, {\Delta }_{12}^{\left(2\right)}=1275с, {\epsilon }_{\min }^{\left(2\right)}=\mathrm{13,76}°C.$ Найденные координаты внутренних точек, в которых достигается минимаксное отклонение, можно видеть на полученной кривой температурного отклонения от заданной величины, представленной на рис. 6.

На рис. 7 представлен график изменения динамики температурного распределения и управляющего воздействия в процессе оптимального по быстродействию индукционного нагрева с ограничением на максимальную температуру точек нагреваемой заготовки. Наглядно виден участок, на котором происходит скачкообразное уменьшение $u_{2\max }$ с целью удержания температуры на уровне $Q_{\text{доп}}$ с того момента ${\upsilon }_0$, когда она впервые достигает допустимого предела при $Q_{\max }>Q_{\text{доп}}$

 

Рис. 7. Изменение управляемой величины и управляющих воздействий с учетом фазового ограничения при коррекции управления ${\mathit{u}}_{\mathbf{2}\mathbf{m}\mathbf{a}\mathbf{x}}$

 

Рассмотрим вариант с ограничением на максимальную мощность первого источника тепла $u_{1\max }$ на интервале интенсивного нагрева неограниченной пластины. Логика построения алгоритма управления не будет иметь серьезных изменений по сравнению с предыдущим случаем. Главное отличие будет заключаться в коэффициенте $L_n$ и в описании конечного температурного поля, то есть вместо (23) и (13) запишем:

$L_n={\theta }_n{A}_n{\eta }_n^2-\frac{u_{2\max }}{u_{1\max }};$ (25)

 

$\begin{array}{c}Q(x,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)})=\sum _{n=1}^{\infty }\frac{\lambda R^2}{a{E}_n{\eta }_n^2}\left(\cos \left({\eta }_n\frac{x}{R}\right)+\frac{B{i}_1}{{\eta }_n}\sin \left({\eta }_n\frac{x}{R}\right)\right)\times \\ \times \left[W_{1n}\left(-1+e^{{\eta }_n^2{\upsilon }_0}+a_Q\left(e^{{\eta }_n^2{\Delta }_{11}^{\left(2\right)}}-e^{{\eta }_n^2{\upsilon }_0}\right)+\frac{b_Q{\eta }_n^2}{{\eta }_n^2-{\beta }_Q}\left(e^{{\eta }_n^2-{\beta }_Q\left({\Delta }_{11}^{\left(2\right)}-{\upsilon }_0\right)}-e^{{\eta }_n^2{\upsilon }_0}\right)\right)\right.\times \\ \times e^{-{\eta }_n^2\left({\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)}\right)}\left.+W_{2n}{u}_{2\max }\left(e^{-{\eta }_n^2\frac{a}{R^2}{\Delta }_{22}^{\left(2\right)})}-e^{-{\eta }_n^2\frac{a}{R^2}({\Delta }_{21}^{\left(2\right)}+{\Delta }_{22}^{\left(2\right)})}\right)\right].\end{array}$ (26)

 

Далее по схеме альтернансного метода остается решить систему (20, б), где $Q\left(x_j^0,{\Delta }_{ij}^{\left(2\right)}\right),i,j=\overline{\mathrm{1,2}}$ определяется уже выражением (26).

Результаты расчета оптимального по быстродействию процесса индукционного нагрева заготовки с ограничением на максимальную температуру при коррекции управления $u_{1\max }$: ${\Delta }_{21}^{\left(2\right)}=1681с, {\Delta }_{22}^{\left(2\right)}=318с, {\Delta }_{11}^{\left(2\right)}=1357с, {\epsilon }_{\min }^{\left(2\right)}=\mathrm{13,35}°C.$Найденные координаты внутренних точек, в которых достигается минимаксное отклонение, можно видеть на полученной кривой температурного отклонения от заданной величины (рис. 8), а графики изменения температуры металлической заготовки в процессе нагрева, а также изменение мощности обоих управляющих воздействий – на рис. 9.

 

Рис. 8. Кривая температурного отклонения при $\mathit{\epsilon }\mathbf{=}{\mathit{\epsilon }}_{\mathbf{m}\mathbf{i}\mathbf{n}}^{\mathbf{\left(}\mathbf{2}\mathbf{\right)}}$, полученная при расчетах с учетом фазового ограничения при коррекции управления ${\mathit{u}}_{\mathbf{1}\mathbf{m}\mathbf{a}\mathbf{x}}$

 

Из двух полученных вариантов решения ЗОУ первый является предпочтительным в силу меньшего времени процесса нагрева.

Существует третий вариант процедуры коррекции режима оптимального управления при поддержании максимальной температуры $Q_{\max }$ на предельно допустимом уровне $Q_{\text{доп}}$ с постепенным уменьшением мощности обоих управляющих воздействий. Каждому из управлений в таком случае должен соответствовать некоторый весовой коэффициент исходя из понимания о соотношении этих составляющих. Это отдельная трудоемкая задача, рассмотрение которой не будет приведено в рамках данной работы.

 

Рис. 9. Изменение управляемой величины и управляющих воздействий с учетом фазового ограничения при коррекции управления ${\mathit{u}}_{\mathbf{1}\mathbf{m}\mathbf{a}\mathbf{x}}$

 

Таким образом, методика расчета оптимального процесса нагрева с учетом технологического ограничения на поведение температурного поля на протяжении всего процесса реализуется в два этапа: на первом этапе необходимо решить задачу быстродействия без учета фазового ограничения и по найденному алгоритму управления построить зависимость $Q(x,{\Delta }_{21}^{\left(2\right)},{\Delta }_{22}^{\left(2\right)},{\Delta }_{11}^{\left(2\right)}),$ по которой проверяется выполнение условий (6). Если условие (6) не выполняется, то алгоритм оптимального по быстродействию управления корректируется по выражениям (21) - (24).

About the authors

Natalya A. Il`ina

Samara State Technical University

Author for correspondence.
Email: ilina.natalyaa@yandex.ru

Graduate student

Russian Federation, 244, Molodogvardeyskaya st., Samara, 443100

References

Supplementary files