Содержание
1. Техническое задание
2. Анализ технического задания
2.1 Прямая задача ВТК
2.2 Обратная задача ВТК
2.3 Модель задачи
2.4 Анализ литературы
2.4.1 Зарубежные методы решения
2.4.2 Отечественные методы решения
3. Прямая задача ВТК для НВТП
3.1 Уравнение Гельмгольца для векторного потенциала
3.2 Поле витка над многослойной средой
3.3 Воздействие проводящего ОК на НВТП
4. Обратная задача ВТК для НВТП
5. Некорректные задачи
5.1 Основные определения. Корректность по Адамару
5.2 Корректность по Тихонову
5.3 Вариационные методы решения некорректных задач
5.3.1 Метод регуляризации
5.3.2 Метод квазирешений
5.3.3 Метод невязки
6. Нелинейное программирование
6.1 Метод штрафных функций
6.2 Релаксационные методы
6.2.1 Метод условного градиента
6.2.2 Метод проекции градиента
6.2.3 Метод случайного спуска
6.3 Метод множителей Лагранжа
7. Линейное программирование
7.1 Алгоритм симплексного метода
8. Одномерная минимизация
8.1 Алгоритм методов
9. Результаты численного моделирования
9.1 Аппроксимации при численном моделировании
9.2 Модели реальных распределений электропроводности
9.3 Принципиальная возможность восстановления
9.4 Восстановление по зашумленным данным
9.5 Восстановление с учетом дополнительной информации
9.6 Восстановление при различном возбуждении
10. Заключение
11. Литература
Приложение 1 — Программная реализация
Приложение 2 — Удельная электропроводность материалов
Приложение 3 — Результаты восстановления
Приложение 4 — Abstract
1. Техническое задание
Разработать алгоритм решения обратной задачи вихретокового контроля (ВТК). Объектом контроля (ОК) являются проводящие немагнитные листы. Объекты контроля подвергаются термообработке (закалка, отпуск) или насыщению внешних слоев различными веществами, что приводит к изменению механических, а вследствие этого и электромагнитных свойств материала листа по глубине.
Задача заключается в определении, в рамках допустимой погрешности, зависимости электропроводности (ЭП) от глубины s (Н) в ОК для данного состояния. Метод контроля заключается в измерении определенного количества комплексных значений вносимой ЭДС на различных частотах с помощью накладного вихретокового преобразователя (НВТП).
Необходимо выбрать математическую модель задачи, способ аппроксимации искомого решения, рассмотреть алгоритм решения.
Используя программную реализацию, исследовать поведение погрешности аппроксимации зависимости s (Н) от следующих факторов:
- От величины приборной погрешности измерения ЭДС
- От вида зависимости электропроводности от глубины s (Н)
- От параметров аппроксимации решения
- От диапазона частот возбуждения ВТП
2. Анализ технического задания.
Основная задача вихретокового контроля с помощью накладных преобразователей состоит из двух подзадач:
- Прямой задачи расчета вносимой ЭДС в присутствии немагнитного проводящего листа с произвольной зависимостью ЭП по глубине.
- Обратной задачи нахождения зависимости ЭП как функции глубины в немагнитном проводящем листе по результатам измерений определенного количества комплексных значений вносимой ЭДС.
2.1 Прямая задача ВТК
Полагая зависимость ЭП от глубины известной проведем ее кусочно-постоянную аппроксимацию. Это позволяет свести исходную задачу к расчету ЭДС в многослойном листе, в каждом слое которого ЭП принимает постоянное значение.
Как показано в работе [50], подобная модель вполне адекватно описывает задачу и дает отличное согласование с результатами опытов.
Рекуррентные формулы для произвольного количества слоев хорошо известны [1-5,36, 42,43,50-52]. Таким образом решение прямой задачи в рамках принятой модели затруднений не вызывает.
2.2 Обратная задача ВТК
С математической точки зрения обратная задача ВТК относится к классу некорректных задач[49] и ее решение неустойчиво т.е. при сколь угодно малой погрешности исходных данных( набора измеренных вносимых ЭДС ) погрешность решения ( рассчитанных локальных значений ЭП ) может быть сколь угодно большой, а одному набору измерений может отвечать много (формально бесконечно много) распределений ЭП по глубине.
При попытке расчета некорректной задачи как корректной, вычислительный процесс за счет неустойчивости сваливается в заведомо худшую сторону. В нашем случае это означает получение распределения ЭП, которое, хотя и обеспечивает требуемое совпадение измеренной и вычисленной ЭДС, но является явно нереальным из-за осцилляций. Следует отметить, что амплитуда и частота осцилляций распределения ЭП растут при увеличении числа независимых параметров аппроксимации ЭП ( коэффициентов полинома в случае полиномиальной аппроксимации, количества узлов при сплайн-аппроксимации и т.д.).
При наличии погрешности измерения вносимой ЭДС, превышающей на несколько порядков вычислительную погрешность и на практике составляющей не менее (0.5-1)% от измеряемого сигнала, ситуация значительно осложняется.
Учитывая вышеизложенное для выделения из множества допустимых распределений решения, наиболее удовлетворяющего физической реальности, в алгоритмах решения обратной задачи необходимо использовать дополнительную априорную информацию. На практике это реализуется введением некоторых критериев, позволяющих отличить решение, отвечающее практике, от физически нереального.
Для решения обратной задачи ВТК предлагались три возможные стратегии[46]:
- Решение большого числа прямых задач и табуляция результатов для различных моделей. Измеренные данные с помощью некоторых критериев сравниваются с таблицей. Подход очень экстенсивный и требующий проведения избыточного числа расчетов, поэтому на практике встречающийся редко.
- Условная минимизация невязки измеренных и расчитанных данных. Очень мощный и универсальный метод, широко распространен для решения обратных задач в различных областях техники [41,44,49]. Позволяет восстанавливать произвольное распределение ЭП по глубине (вообще говоря произвольное 3D распределение), но требуется довольно сложная процедура расчета.
- Аналитическое инвертирование ядра оператора и использование алгоритма, зависящего от ядра уравнения[46]. Потенциально самый малозатратный метод, однако как и все аналитические, применим далеко не всегда.
В нашем случае остановимся на втором подходе, поскольку он сочетает в себе универсальность, точность и относительную простоту реализации.
В целом процесс решения обратной задачи сводится к итерационному решению прямой задачи для текущей оценки распределения ЭП и внесению изменений в эту оценку в соответствии с величиной невязки.
2.3 Модель задачи
Приведем основные положения, на основе которых будет построена модель нашей задачи:
- ОК представляет из себя находящуюся в воздухе проводящую пластину толщиной Н состоящую из N плоско-параллельных слоев толщиной b i .
- В пределах каждого слоя удельная электропроводность s имеет постоянное значение т.е. распределение s по глубине аппроксимируется кусочно-постоянной зависимостью.
- Возбуждающая и измерительная обмотки ВТП заменяются нитевидными моделями. Следует отметить, что это предположение сказывается лишь на решении прямой задачи, а проведя интегрирование можно получить выражения для катушек конечных размеров.
- Для численного моделирования реальных распределений ЭП применим пять типов аппроксимации: сплайном, кусочно-постоянную, кусочно-линейную, экспоненциальную и гиперболическим тангенсом. В процессе решения прямой задачи с их помощью вычисляются значения s в центральных точках слоев пластины.
2.4 Анализ литературы
2.4.1 Зарубежные методы решения
Решению обратной задачи ВТК посвящен ряд работ в зарубежных изданиях. Следует отметить монографию [38], в которой рассмотрены случаи импульсного возбуждения, а оперируют в частотной и временной областях напряженностью электрического поля.
Подход к решению квазистационарных задач рассмотрен в цикле статей [45-51]. Он основан на интегральной постановке задачи с помощью функций Грина[31-34,39]. Для иллюстрации рассмотрим решение обратной задачи ВТК согласно [49].
А. Прямая задача
Определим функцию v(r)=( s (r) — s 0 )/ s 0 , где s (r) — произвольное распределение проводимости, а s 0 — ее базовая величина. Функция v(r) может представлять собой как описание произвольного распределения проводимости (в этом случае для удобства полагаем s (r)= s 0 вне некоторого ОК объема V , тогда v(r) отлична от нуля только в пределах V ) так и некоторого дефекта (для трещины v(r)=-1 внутри дефекта и равна нулю вне его).
Рассмотрим систему уравнений Максвелла в предположении гармонического возбуждения exp(-jwt) и пренебрегая токами смещения:
|
|
( 2.4.1) |
где P(r)=[ s (r)- s 0 ] Ч E(r)= s 0 Ч v(r) Ч E(r) — может интерпретироваться как плотность диполей эффективного тока, причиной которого является вариация s (r)- s 0 .
Решение уравнений Максвелла можно представить в виде
|
|
( 2.4.2) |
где E i (r) — возбуждающее поле, а G(r|r’) — функция Грина, удовлетворяющая уравнению С ґ С ґ G(r|r’)+k 2 Ч G(r|r’)= d (r-r’) , k 2 =-j Ч w Ч m 0 Ч s 0 , d (r-r’) — трехмерная дельта-функция.
Импеданс ВТП можно выразить как
|
|
( 2.4.3) |
где интеграл берется по измерительной катушке, J(r) — плотность тока в возбуждающей катушке. Применяя теорему взаимности импеданс можно представить через возбуждающее поле:
|
|
( 2.4.4) |
где интеграл берется по объему ОК.
В. Обратная задача
Пусть v(r) — оценка истинной функции v true (r) , Z obs (m) — измеренный импеданс ВТП в точке r 0 на частоте возбуждения w , m=(r 0 , w ) — вектор в некоторой области определения M , Z[m,v] — оценка величины Z obs (m) на основе решения прямой задачи.
Определим функционал невязки измеренных и рассчитанных значений импеданса ВТП как :
|
|
( 2.4.5) |
Предположим, что для решения обратной задачи используется итерационный алгоритм типа метода спуска: v n (r)= v n-1 (r)+ a s n (r) . Можно показать, что в случае метода наискорейшего спуска итерация имеет вид: v n (r)= v n-1 (r)- a Ч С F[ v n-1 (r) ] , где градиент функционала С F[v] можно определить как :
|
|
( 2.4.6) |
где Re обозначает вещественную часть, * обозначает комплексную сопряженность.
Требуемый в (2.4.6) градиент импеданса можно определить как:
|
С Z(r) = — s 0 Ч E(r) Ч E * (r) |
( 2.4.7) |
где E * (r) — решение уравнения
|
|
( 2.4.8) |
С. Аппроксимация при решении обратной задачи
Пусть электропроводность моделируется с помощью конечного числа переменных (например узловых значений некоторой аппроксимации), а вектор р состоит из этих переменных. Тогда выражение (2.4.7) принимает вид:
|
|
( 2.4.9) |
где ( С Z) j — j -ая компонента градиента импеданса.
Значение j-ой компоненты градиента невязки (2.4.6) можно представить как:
|
|
( 2.4.10) |
Следует обратить внимание на то, что в случае дискретного пространства М (конечное число измерений) интеграл в (2.4.10) заменяется суммой.
С учетом приведенных преобразований итерация метода наискорейшего спуска принимает вид:
|
p j n = p j n-1 — a Ч ( С F n-1 ) j |
( 2.4.11) |
где n — номер итерации.
D. Пример применения
В качестве примера рассмотрим функцию v(r) в виде v(r)= S c i Ч f i (r) , i=1,N , где f i (r) — множество линейно независимых базовых функций с коэффициентами c i . Рассматривая коэффициенты c i в роли параметров аппроксимации ( c i =p i ) получим из (2.4.9) для компонентов градиента импеданса:
|
|
( 2.4.12) |
В случае проводящего ОК, состоящего из N параллельных слоев с проводимостью s j распределение электропроводности по глубине можно представить с помощью функций Хевисайда H(z) как s (z)= S s j Ч [ H( z-z j ) — H( z-z j+1 ) ] .
Подставляя в (2.4.12) базовые функции вида f i (z)=[H( z-z j )-H( z-z j+1 )] , получим окончательное выражение:
|
|
( 2.4.13) |
Отметим основное преимущество такого решения. Несмотря на определенную сложность вычислений при решении интегральных уравнений (2.4.2-2.4.8) для расчета градиента импеданса НВТП необходимо решить только две такие задачи.
2.4.2 Отечественные методы решения
Подход, в значительной мере аналогичный работам [45-51] был предложен в работе [41]. Из-за небольшого объема в ней уделено недостсточное внимание вопросам практической реализации, объяснены не все обозначения и не приведены результаты численного моделирования. В целом это значительно снижает практическую ценность статьи. Приведем основные положения этой работы.
Прямая задача
Пусть круговой виток радиусом а с током I находится в точке P=P s (r, j ,z), j О (- p , p ) вблизи немагнитного ОК, занимающего область V . Пусть ОК обладает электрической проводимостью s = s 0 Ч s (Р) являющейся произвольной функцией координат. Требуется по N измерениям величины э.д.с. определить s как функцию координат точек P О V . Причем i -ое измерение э.д.с. будем проводить на i -ом измерительном круговом витке с координатами P i =P i (r, j ,z) i=1,N при неизменных частоте и расположении возбуждающего витка.
В общем случае напряженность электрического поля Е определяется через векторный магнитный потенциал А , причем А = А 0 + А вн , где А 0 — возбуждающий, а А вн — вносимый потенциалы.
|
|
(2.4.14) |
Вводя функцию Грина G(p,p 0 ) получим
|
|
(2.4.15) |
При этом вносимая напряженность электрического поля
|
E вн = -j Ч w Ч A вн |
(2.4.16) |
Вносимая э.д.с., наводимая в i -ом витке
|
|
(2.4.17) |
где функция Грина G(P,P 0 ) имеет вид
|
|
(2.4.18) |
В дальнейшем рассмотрим случай, при котором V -полупространство (r>0,\ j \< p ,z 0
(5.2)
Используя классический регуляризирующий функционал вида 
в терминах нашей задачи получаем:
|
|
(5.3) |
Основное преимущество метода состоит в регуляризации простейшим способом, в рамках использования квадратичного функционала. Это позволяет использовать для решения некорректной задачи хорошо известные и легко программируемые методы минимизации квадратичных функционалов [17].
Оборотной стороной достоинств метода являются его недостатки. Требование минимизации нормы решения и, как следствие, выбор гладкой реализации, в нашем случае будет приходить в противоречие с физикой задачи и в принципе не позволит находить решения с выраженным приповерхностным изменением ЭП.
Еще один принципиальный недостаток метода состоит в постановке функционала как квадратичного, единого для всех измерений. Его минимум в общем случае не гарантирует минимизацию отклонения для произвольного i -го измерения в следствии нелокальности условия минимизации.
Кроме того, следует учитывать отсутствие надежных априорных рекомендаций по выбору параметра регуляризации a . Обычно подходящие значения a можно подобрать только после ряда численных экспериментов по решению однотипных задач. Изменение характера искомого решения приводит к необходимости поиска нового значения a .
5.3.2 Метод квазирешений
Метод использует одну из форм критерия невязки и заключается в сведении невязки к минимуму на некотором непустом множестве P , содержащем подмножество искомых решений.
Квазирешением уравнения (5.1) на множестве P М X называется всякий элемент y О P для которого справедливо равенство r F ( Ay , f ) = inf( Ax , f ), x О P . Понятие квазирешения обобщает понятие решения, а для его существования не требуется принадлежность решения множеству P .
Исходя из вышеизложенного получаем постановку метода в виде задачи условной минимизации функционала Ф(x,f) :
|
|
(5.4) |
Отметим, что множество Р может иметь простой вид, например интервала [ x min , x max ].
В терминах нашей задачи ВТК постановка задачи (5.4) примет вид:
|
|
(5.5) |
Для того, чтобы гарантировать минимизацию отклонения для произвольного i -го измерения, можно применить к первому выражению в (5.5) локальный в смысле Чебышева критерий, в соответствии с которым получаем окончательное выражение :
|
|
(5.6) |
Основное преимущество метода состоит в том, что само понятие квазирешения снимает трудности с требованиями тихоновской корректности: первым (вызывающим переопределенность задачи) и третьим (обычно принадлежность приближенной правой части уравнения (5.1) множеству N=AM неизвестна, а критерии этой принадлежности часто сами бывают неустойчивы).
Кроме этого при рассмотрении задачи в виде (5.6) возможна постановка минимизационной задачи как задачи нелинейного программирования с явно заданными ограничениями на искомые переменные. В этом случае нет необходимости искажать исходный функционал регуляризующими членами как в (п5.3.1), а требования к искомому решению можно удовлетворить, управляя ограничениями на параметры минимизации (в нашем случае — узловые значения ЭП).
5.3.3 Метод невязки
Рассмотрим множество Р формальных решений уравнения (5.1) Р={x : r F ( Ax , f d ) Ј d }, где f d — приближенная правая часть (5.1), известная с погрешностью d .
В качестве приближенного решения (5.1) нельзя брать произвольный элемент множества Р , т.к. не гарантируется близость Р к множеству точных решений. Для выбора приближенного решения предлагается использовать стабилизирующий функционал W (х) из (п 5.3.1) следующим образом: W ( х ) = inf W ( х ), x О P . Этот способ приводит к выбору элементов множества Р имеющих минимально допустимую невязку.
С учетом этого постановка метода состоит в условной минимизации функционала Ф(х) :
|
|
(5.7) |
Как и для метода регуляризации можно использовать стабилизирующий функционал вида W (х)=||x|| 2 , что приводит в обозначениях нашей задачи к системе:
|
|
(5.8) |
При использовании локального в смысле Чебышева критерия система (5.8) окончательно примет вид:
|
|
(5.9) |
6. Нелинейное программирование
Содержание нелинейного программирования составляют теория и методы решения задач о нахождении экстремумов нелинейных функций на множествах, определяемых линейными и нелинейными ограничениями (равенствами и неравенствами)[16-29].
Рассмотрим наиболее распространенные методы решения на примере основной задачи нелинейного программирования вида:
|
|
(6.1) |
6.1 Метод штрафных функций
Идея метода состоит в замене экстремальной задачи с ограничениями (6.1) на задачу безусловной минимизации однопараметрической функции
|
|
(6.2) |
, b >0Непрерывную функцию y (х) называют штрафом, если y (х)=0 для х О Х и y (х)>0 в противном случае. Функция y (х) должна быть выбрана таким образом, чтобы решение задачи (6.2) сходилось при b ® 0 к решению исходной задачи (6.1) или, по крайней мере, стремилось к нему.
Приведем часто используемые выражения для штрафа :
|
|
(6.3) |
|
|
(6.4) |
|
|
(6.5) |
, k>0
Наибольшее применение находит штраф (6.3). Выражение (6.5) гарантирует конечность метода при любом k>0 .
При численной реализации метода штрафных функций возникают проблемы выбора начального значения параметра b и способа его изменения. Сложность состоит в том, что выбор достаточно малого b увеличивает вероятность сходимости решения (6.2) к решению (6.1), а скорость сходимости градиентных методов вычисления точек минимума (6.2), как правило, падает с убыванием величины b .
6.2 Релаксационные методы
Релаксационным методом называют процесс построения последовательности точек {х k : х k О X , j ( х k+1 ) Ј j ( х k ) ; k=0,1… }. Основными представителями этого класса являются методы спуска, алгоритм которых состоит из следующих шагов :
- Выбор начального приближения х 0
- Выбор в точке х k направления спуска -s k
- Нахождение очередного приближения х k+1 = х k — a k Ч s k , где длина шага a k >0
Различия методов состоят в выборе либо направления спуска, либо способа движения вдоль выбранного направления. В последнем случае обычно используют одномерную минимизацию функции х k+1 ( a ) = х k — a Ч s k (при этом точность вычисления точки минимума функции х k+1 ( a ) следует согласовывать с точностью вычисления значений функции j (х) ) или способ удвоения a (величина шага удваивается пока выполняется условие j (х k+1 ) Ј j (х k ) ).
6.2.1 Метод условного градиента
Идея метода заключается в линеаризации нелинейной функции j (х). В этом методе выбор направления спуска осуществляется следующим образом :
- Линеаризируя функцию j (х) в точке х К получаем Ф(х)= j ( х k ) + ( j ( х k )’ , х — х k )
- Минимизируя линейную функцию Ф(х) на множестве Х находим х min
- Направление спуска получаем как -s k = х min — х k
Таким образом итерация метода имеет вид: x k+1 =x k + a k Ч (s k+1 — x k ) , s k+1 =arg min( С f(x k ),x).
Основное преимущество метода проявляется в случае задания допустимого множества с помощью линейных ограничений. В этом случае получаем задачу линейного программирования, решаемую стандартными методами(например симплексным).
6.2.2 Метод проекции градиента
Этот метод является аналогом метода градиентного спуска, используемого в задачах без ограничений. Его идея состоит в проектировании точек, найденных методом наискорейшего спуска, на допустимое множество, определяемое ограничениями. Проекцией точки y на множество Х называется точка P(y) О Х такая, что || P(y) — y || Ј || x — y || для всех х О Х. Задача проектирования формализуется как || x — y || 2 ® min, x О Х.
Выбор направления спуска осуществляется следующим образом :
- Находим точку r k = х k — a Ч j ’( х k )
- Находим проекцию p k точки r k на множество Х
- Направление спуска получаем как -s k = p k — х k
Таким образом итерация метода имеет вид: x k+1 =P X [ x k — a k Ч С f( x k ) ], где Р X (у) — ортогональная проекция точки у на множество Х .
Для отыскания направления спуска s k необходимо решить задачу минимизации квадратичной функции || r k — х || 2 на множестве Х . В общем случае эта задача того же порядка сложности, что и исходная, однако для задач, допустимое множество которых имеет простую геометрическую структуру, отыскание проекции значительно упрщается. Например, для многомерного параллелепипида Q N ={x О R N : a Ј x Ј b }, отыскание проекции осуществляется путем сравнения n чисел и имеет вид P(x)={ a i , x i b i }.
6.2.3 Метод случайного спуска
Метод характеризуется тем, что в качестве направления спуска s K выбирается некоторая реализация n -мерной случайной величины S с известным законом распределения. Об эффективности этого метода судить трудно, однако благодаря использованию быстродействующих ЭВМ он оказывается практически полезным.
6.3 Метод множителей Лагранжа
Идея метода состоит в отыскании седловой точки функции Лагранжа задачи (6.1). Для нахождения решения вводится набор переменных l i , называемых множители Лагранжа, и составляется функция Лагранжа, имеющая вид:
|
|
(6.6) |
Алгоритм метода состоит в следующем:
- Составление функции Лагранжа
- Нахождение частных производных функции Лагранжа
|
|
(6.7) |
- Решение системы из n+m уравнений вида
|
|
(6.8) |
Решениями системы (6.8) являются точки, которые могут быть решениями задачи.
- Выбор точек, в которых достигается экстремум и вычисление функции j (х) в этих точках.
7. Линейное программирование
Задача линейного программирования в каноническом виде имеет вид[15,16]:
|
|
(7.1) |
Приведение к каноническому виду любой задачи линейного программирования осуществляется путем введения дополнительных неотрицательных переменных, за счет чего ограничения, имеющие вид неравенств, принимают вид эквивалентных им равенств.
Любая задача линейного программирования может быть решена за конечное число итераций с помощью симплексного метода[17,18]. Следует отметить, что поскольку этот метод разработан для неотрицательных элементов x j , это условие учитывается неявно и в систему уравнений (7.1) при численной реализации не входит.
7.1 Алгоритм симплексного метода
1. Приведение к каноническому виду
2. Выбор начального базиса
3. Проверка оптимальности базиса
Матрицу А можно рассматривать как совокупность столбцов a j т.е. е a j Ч x j =b где j=1,N . Не ограничивая общности можно считать, что базис образуют первые m столбцов, тогда остальные можно представить в виде a k = е a j Ч l jk , j=1,m где l jk . — некоторые числа.
Рассмотрим коэффициенты D k = е c j Ч l jk — c k где j=1,m и k=1,N . Заметим, что для базовых столбцов D k є 0 . Проверка на оптимальность осуществляется следующим образом:
|
D k < 0 , k=1,N |
— текущий базис оптимален |
|
|
— решение не ограничено сверху |
|
|
— существует другой, более подходящий базис |
4. Составление нового базиса
- Если ( h k Ј e ) то x min =min{ j (c k ) , j (d k ) } иначе k++ и переход к шагу II
4.1 Выбор элемента для введения в базис.
В базис вводится любой столбец, для которого D k j ( d k )
a k =
a k-1
c k-1
b k =
d k-1
b k-1
d k =
c k-1
c k =
d k-1
h k+2 =
h k -h k-1
h k -h k-1
d k =
b k -h k+2
c k =
a k +h k+2
Следует отметить, что на каждом шаге кроме первого, производится только одно вычисление значения функции j (x) .
Легко показать, что для получения оптимальной последовательности отрезков, стягивающихся к точке минимума, необходимо положить l k = F k-1 / F k , где F — число Фибоначчи.
8.2 Метод Фибоначчи
Решая вопрос, при каких значениях параметра l за конечное число итераций N мы получим отрезок минимальной длины, получим l = l N = F N-1 / F N . Иначе говоря, для поиска минимума первоначально необходимо найти число Фибоначчи N такое, что F N+1 0000XXXX}
fHypTg:=(( apDT AND apHypTg ) = apHypTg);
if fHypTg then
begin
si0[ 1 ]:=si[ 1 ]; {si1 — conductivity about bottom of slab}
si0[ 2 ]:=par0[ 2 ]; {si2 — conductivity about top of slab}
si0[ 3 ]:=par0[ 3 ]; {Beta — ratio of approx.}
si0[ 4 ]:=par0[ 4 ]; {Gamma- ratio of approx.}
mCur:=4;
end
else
if(( apDT AND apExp ) = 0 ) then {It’s not an EXP approx.}
begin
for i:=1 to nPoints do si0[ i ] :=si [ i ]; {SI data from file}
mCur:=nPoints;
end
else
begin
si0[ 1 ]:=si[ 1 ]; {siI — conductivity about bottom of slab}
si0[ 2 ]:=si[ nPoints ]; {siE — conductivity about top of slab}
si0[ 3 ]:=par0[ 1 ]; {Alfa- ratio of approx.}
mCur:=3;
end;
setApproximationType( apDT ); {approx. type for direct problem}
setApproximationData( si0, mCur ); {approx. data for direct problem}
nApprox := ( nApprox AND $0F ); {XXXXYYYY->0000YYYY}
fHypTg := (( nApprox AND apHypTg ) = apHypTg );
fMulti := (( nApprox AND apExp ) = 0 ) AND NOT fHypTg; {It’s not an EXP approx.}
if fMulti then
begin
for i:=1 to nPoints do
begin
Gr[ 1,i ]:=SiMax[ i ];
Gr[ 2,i ]:=SiMin[ i ];
Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2; {zero estimate of SI}
Rgs[ i ]:=1E33; {biggest integer}
end;
mLast:=nPoints; {loop for every node of approx.}
mCur :=1; {to begin from the only node of approx}
end
else
if fHypTg then
begin
Gr[ 1,1 ]:= siMax[ 1 ]; Gr[ 2,1 ]:= siMin[ 1 ]; Rgs[ 1 ]:=1E33;
Gr[ 1,2 ]:=parMax[ 2 ]; Gr[ 2,2 ]:=parMin[ 2 ]; Rgs[ 2 ]:=1E33;
Gr[ 1,3 ]:=parMax[ 3 ]; Gr[ 2,3 ]:=parMin[ 3 ]; Rgs[ 3 ]:=1E33;
Gr[ 1,4 ]:=parMax[ 4 ]; Gr[ 2,4 ]:=parMin[ 4 ]; Rgs[ 4 ]:=1E33;
for i:=1 to 4 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;
mLast:=1;
mCur:=4;
end
else
begin
Gr[ 1,1 ]:= siMax[1]; Gr[2,1]:= siMin[1]; Rgs[ 1 ]:=1E33;
Gr[ 1,2 ]:= siMax[nPoints]; Gr[2,2]:= siMin[nPoints]; Rgs[ 2 ]:=1E33;
Gr[ 1,3 ]:= parMax[1]; Gr[2,3]:= parMin[1]; Rgs[ 3 ]:=1E33;
for i:=1 to 3 do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;
mLast:=1;
mCur :=3;
end;
initConst( nLayers, parMaxH, parMaxX , parEps, parEqlB );{set probe params}
end;
procedure directTask; {emulate voltage measurements [with error]}
begin
for i:=1 to nFreqs do
begin
getVoltage( freqs[i], Umr[ i ], Umi[ i ] ); {«measured» Uvn*}
if ( epsU > 0 ) then {add measurement error}
begin
randomize; Umr[ i ]:=Umr[ i ]*( 1 + epsU*( random-0.5 ) );
randomize; Umi[ i ]:=Umi[ i ]*( 1 + epsU*( random-0.5 ) );
end;
end;
writeln(‘* Voltage measurements have been emulated’);
setApproximationType( nApprox ); {approx. type for inverse problem}
setApproximationData( Rg, mCur ); {approx. data for inverse problem}
end;
procedure reduceSILimits; {evaluate SI for m+1 points of approx. using aG}
var
x0, x1, xL, dx, Gr1, Gr2 : real;
j, k : byte;
begin
{—————————— get SI min/max for m+1 points of approximation}
dx:=1/( nPoints-1 );
for i:=1 to m+1 do
begin
k:=1;
x1:=0;
x0:=( i-1 )/m;
for j:=1 to nPoints-1 do
begin
xL:=( j-1 )/( nPoints-1 );
if( ( xL < x0 ) AND ( x0 Gr[1,i] )then Rg[i]:=Gr[1,i];
if ( Rg[i] < Gr[2,i] )then Rg[i]:=Gr[2,i];
if m > 1 then {There’re more than 1 point of approx.}
begin
Gr1:= Rg[i]+( Gr[1,i]-Rg[i] )*aG; {reduce upper bound}
Gr2:= Rg[i]-( Rg[i]-Gr[2,i] )*aG; {reduce lower bound}
if ( Gr1 < Gr[1,i] )then Gr[1,i]:=Gr1; {test overflow}
if ( Gr2 > Gr[2,i] )then Gr[2,i]:=Gr2;
end;
end;
setApproximationData( Rg , m+1 );
end;
procedure resultMessage; {to announce new results}
begin
if fMulti then
begin
writeln(‘ current nodal values of conductivity’);
write(‘ si : ‘);for i:=1 to m do write(Rg[i] :6:3,’ ‘);writeln;
write(‘ max: ‘);for i:=1 to m do write(Gr[1,i]:6:3,’ ‘);writeln;
write(‘ min: ‘);for i:=1 to m do write(Gr[2,i]:6:3,’ ‘);writeln;
end
else
begin
for i:=1 to nPoints do si[i]:=getSiFunction( ( i-1 )/( nPoints-1 ) );
if fHypTg then
saveHypTgResults
else
saveExpResults;
end;
end;
procedure clockMessage; {user-friendly message}
begin
writeln(‘***********************************************************’);
write( ‘* approximation points number :’,m:3,’ * Time ‘); clock;
writeln(‘***********************************************************’);
end;
procedure done; {final message}
begin
Sound(222); Delay(111); Sound(444); Delay(111); NoSound; {beep}
write(‘* Task processing time ‘); clock; saveTime;
writeln(‘* Status: Inverse problem has been successfully evaluated.’);
end;
Begin
about;
loadData;
initParameters;
directTask;
for m:=1 to mLast do
begin
if fMulti then
begin
mCur:=m;
clockMessage;
end;
doMinimization; {main part of work}
setApproximationData( Rg, mCur ); {set new approx. data}
resultMessage;
if(( fMulti )AND( m < nPoints ))then reduceSILimits;
end;
done;
End.
П1.5 Модуль глобальных данных EData
Unit EData;
Interface
Uses DOS;
Const
maxPAR = 40; {nodes of approximation max number}
maxFUN = 20; {excitation frequencies max number}
maxSPC = 4; {support approximation values number}
iterImax = 50; {max number of internal iterations}
Const
apSpline = 1; {approximation type identifiers}
apHypTg = 3;
apExp = 2;
apPWCon = 4;
apPWLin = 8;
Type
Parameters = array[ 1..maxPAR ] of real; {Si,Mu data}
Functionals = array[ 1..maxFUN ] of real; {Voltage data}
SpecialPar = array[ 1..maxSPC ] of real; {Special data}
Var
hThick : real; {thickness of slab}
nPoints : integer; {nodes of approximation number within [ 0,H ]}
nLayers : integer; {number of piecewise constant SI layers within[ 0,H ]}
nFreqs : integer; {number of excitation frequencies}
nStab : integer; {required number of true digits in SI estimation}
epsU : real; {relative error of measurement Uvn*}
nApprox : byte; {approximation type identifier}
incVal : real; {step for numerical differ.}
parMaxH : integer; {max number of integration steps}
parMaxX : real; {upper bound of integration}
parEps : real; {error of integration}
derivType: byte; {1 for right; 2 for central}
Var
freqs : Functionals; {frequencies of excitment Uvn*}
Umr, Umi : Functionals; { Re(Uvn*),Im(Uvn*) for «measured» data}
Uer, Uei : Functionals; { Re(Uvn*),Im(Uvn*) for estimated data}
mu : Parameters; {relative permeability nodal values}
si, si0 : Parameters; {conductivity approximation nodal values}
siMin, siMax : Parameters; {conductivity nodal values restrictions}
par0 : SpecialPar; {alfa,si2,beta,gamma — for exp&HypTg}
parMin,parMax: SpecialPar; {-||- min/max}
zLayer : Parameters; {relative borders of slab layers [0,1]}
Var
aG : real; {scale factor for SImin/max}
Ft : real; {current discrepancy functional value}
fMulti : boolean; {TRUE if it isn’t an EXP-approximation}
fHypTg : boolean; {TRUE for Hyperbolic tg approximation}
parEqlB : boolean; {TRUE if b[i]=const}
mCur : integer; {current number of approximation nodes}
inFileName : string; {data file name}
outFileName : string; {results file name}
Var
Rg : Parameters; {current SI estimation}
RgS : Parameters; {previous SI estimation}
Gr : array [ 1..2 ,1..maxPAR ] of real; {SI max/min}
Fh : array [ 1..maxPAR , 1..maxFUN ] of real; {current discrepancies}
Type
TTime=record
H, M, S, S100 : word; {hour,min,sec,sec/100}
end;
Var
clk1, clk2 : TTime; {start&finish time}
procedure loadData; {load all user-defined data from file}
Implementation
procedure loadData;
var
i,eqlB : integer;
FF : text;
begin
assign( FF, outFileName ); {clear output file}
rewrite( FF );
close( FF );
assign( FF, inFileName ); {read input file}
reset( FF );
readln( FF );
readln( FF );
readln( FF, hThick, nPoints, nLayers, nFreqs, nStab, epsU, aG, nApprox );
readln( FF );
for i:=1 to nFreqs do read( FF, freqs[i] );
readln( FF );
readln( FF );
readln( FF );
for i:=1 to nPoints do readln( FF, si[i], siMin[i], siMax[i] );
readln( FF );
readln( FF );
readln( FF , incVal, parMaxH, parMaxX, parEps, derivType, eqlB );
readln( FF );
readln( FF );
for i:=1 to maxSPC do readln( FF, par0[i] , parMin[i] , parMax[i] );
readln( FF );
if ( eqlB=0 )then
begin
for i:=1 to nLayers+1 do read( FF, zLayer[i] );
parEqlB:=false;
end
else parEqlB:=true;
close( FF );
for i:=1 to maxPAR do mu[i]:=1;
end;
Var
str : string;
Begin
if( ParamCount = 1 )then str:=ParamStr(1)
else
begin
write(‘Enter I/O file name, please: ‘);
readln( str );
end;
inFileName :=str+’.txt’;
outFileName:=str+’.lst’;
End.
П1.6 Модуль работы с файлами EFile
Unit EFile;
Interface
Uses
DOS, EData;
function isStable( ns : integer; var RG1,RG2 ) : boolean;
function saveResults( ns,iter : integer ) : boolean;
procedure saveExpResults;
procedure saveHypTgResults;
procedure clock;
procedure saveTime;
Implementation
Var
FF : text;
i : byte;
function decimalDegree( n:integer ) : real;{10^n}
var
s:real; i:byte;
begin
s:=1;
for i:=1 to n do s:=s*10;
decimalDegree:=s;
end;
function isStable( ns:integer ; var RG1,RG2 ) : boolean;
var
m : real;
R1 : Parameters absolute RG1;
R2 : Parameters absolute RG2;
begin
isStable:=TRUE;
m:=decimalDegree( ns-1 );
for i:=1 to mCur do
begin
if NOT(( ABS( R2[i]-R1[i] )*m ) dx ) do
begin
i:=i + 1;
dx:=dx + dh;
end;
siPWConst:=appSigma[ i ];
end;
end;
function siPWLinear( x:real ) : real;{Piecewise linear approximation}
var
dx, dh : real;
i : byte;
begin
if( appCount = 1 )then siPWLinear := appSigma[ 1 ]
else
begin
dh:=1/( appCount-1 );
dx:=0;
i:=1;
repeat
i:=i + 1;
dx:=dx + dh;
until( x B[k] then k1:=k;
for k:=1 to mp1 do A[k1,k]:=-A[k1,k];
A[k1,mCur+1+k1]:=0;
B[k1]:=-B[k1];
for i:=1 to n2 do
if ik1 then
begin
B[i]:=B[i]+B[k1];
for k:=1 to mm1 do A[i,k]:=A[i,k]+A[k1,k];
end;
for i:=mp2 to m21 do
begin
Sx[i]:=B[i-mp1];
Nb[i-mp1]:=i;
end;
for i:=1 to mp1 do Sx[i]:=0;
Sx[1]:=B[k1];
Sx[mp1+k1]:=0;
Nb[k1]:=1;
103:
for i:=2 to m21 do N0[i]:=0;
104:
for i:=m21 downto 2 do
if N0[i]=0 then n11:=i;
for k:=2 to m21 do
if ((A[k1,n11]B[i]/A[i,n11] then
begin
Sx[n11]:=B[i]/A[i,n11]; ip:=i;
end;
end;
end
else
if iq=0 then
begin
N0[n11]:=n11;
goto 104;
end;
end;
Sx[Nb[ip]]:=0;
Nb[ip]:=n11;
B[ip]:=B[ip]/A[ip,n11];
apn:=A[ip,n11];
for k:=2 to m21 do A[ip,k]:=A[ip,k]/apn;
for i:=1 to m1 do
if iip then
begin
ain:=A[i,n11];
B[i]:=-B[ip]*ain+B[i];
for j:=1 to m21 do A[i,j]:=-ain*A[ip,j]+A[i,j];
end;
for i:=1 to m1 do Sx[Nb[i]]:=B[i];
goto 103;
105:
for k:=1 to mCur do Sx[k+1]:=Sx[k+1]+Gr[2,k];
a1:=0;
a2:=1.;
dh:=a2-a1;
r:=0.618033;
tl:=a1+r*r*dh;
tp:=a1+r*dh;
j:=1;
108:
if j=1 then tt:=tl else tt:=tp;
106:
for i:=1 to mCur do Rg[i]:=Zt[i]+tt*(Sx[i+1]-Zt[i]);
getFunctional( 0 );
cv:=abs(Fh[1,1]);
if nFreqs>1 then
for k:=2 to nFreqs do
begin
cv1:=abs(Fh[1,k]);
if cv