Форум АСУ в Україні

Here a finding of Dahlquist's gives useful guidance. Dahlquist defines an integration algorithm as having the highly desirable property of 'A-stability' if it is stable for all step lengths when applied to the linear differential equation describing an unconditionally stable physical system

Ось рішення Далквіста дає корисні рекомендації. Далквіст визначає алгоритм інтегрування які бажано мають високу властивістю 'A-стабільності ", якщо вона стабільна у всіх довжинах кроку при нанесенні на лінійне диференціальне рівняння описує безумовно стійку фізичну систему

with λ strictly positive. Only implicit algorithms of order two or below are A-stable.

з λ строго позитивні. Тільки неявні алгоритми другого порядку або нижче, А-стабільні.

To illustrate the difference in stability properties between explicit and implicit integration algorithms, consider again the equation used to describe valve dynamics in Section 2.2. Dropping the subscripts from equation (2.9) for clarity and generality, and setting the demanded valve travel, Xd, to zero, indicating a demand for closure, we have:

Щоб проілюструвати різницю у властивостях стабільності у відносинах між явними і неявними алгоритмами інтегрування, знову розглянемо рівняння, яке використовується для опису динаміки клапанів в розділі 2.2. Опустивши індекси з рівняння (2.9) для ясності і узагальнення, а також створення необхідного ходу штоку клапана, Xd, до нуля, що вказує на потребу у закритті, ми повинні:

which is in the form of equation (2.63). The explicit, forward Euler approximation yields the equation

який знаходиться у вигляді рівняння (2,63). Явно, попереднє наближення Ейлера дає рівняння

which leads to the condition stated in the first paragraph of this section, namely

що призводить до умов, зазначених у першому абзаці цього розділу, а саме

However, the implicit, backward Euler approximation gives

Тим не менш, неявна, тому наближення Ейлера дає


From inspection, this is stable for all positive values of the timestep, At, in-line with Dahlquist.

З огляду, це стабільна, для всіх позитивних значень тимчасового кроку, Δt, на лінії Дальквіста.

But while the implicit integration algorithm above can be rearranged easily into an explicit form for the single-variable, linear case, the same cannot be said for the multivariable, nonlinear cases that we will normally be dealing with in process modelling. If we examine the general simulation case given by equation set (2.22), then applying the implicit, backward Euler algorithm produces the set of equations:

Але в той час як неявний алгоритм інтегрування поданий вище, може бути легко переобладнаний в явному вигляді для однієї змінної в лінійному випадку, того ж не можна сказати про багато змінних в нелінійних випадках, з якими ми, як правило, матимемо справу в моделюванні процесів. Якщо ми розглянемо загальний випадок моделювання, який визначається рівнянням в наборі (2,22), то, застосовуючи неявний зворотній алгоритм Ейлера проводить до набіру рівнянь:

which represents a system of n nonlinear simultaneous equations in the n unknowns of the vector x at the (k + 1 )th timestep. We will need to solve a similar set of simultaneous equations at each timestep. Thus in order to get the boon of an algorithm with much better stiffness properties, allowing us to take much bigger timesteps, we have had to pay for it by involving ourselves in significantly more computation at each timestep.

який являє собою систему n нелінійних сумісних рівнянь з n невідомими векторами х на (k + 1) -го тимчасового кроку. Нам потрібно буде вирішити подібну серію сумісних рівнянь на кожному часовому кроці. Таким чином, для того, щоб отримати благо алгоритму з набагато кращими властивостями жорсткості, що дозволяє нам приймати набагато більше тимчасових кроків, ми повинні заплатити за нього, залучаючи себе в значно більшу кількість обчислень на кожному часовому кроці.

We will now illustrate the way that equation (2.70) could be solved as part of a stiff integration package. The solution relies partly on using the Newton-Raphson technique for solving nonlinear simultaneous equations, the principles of which will now be explained. We may describe a system of n nonlinear, simultaneous equations in the n unknowns of the vector z by the vector equation:

Тепер проілюструюємо спосіб, в якому рівняння (2.70) може бути вирішене в рамках жорсткого інтегрування. Рішення спирається частково на використанні методів Ньютона-Рафсона для вирішення нелінійних систем рівнянь, принципи, які тепер будуть пояснені. Ми можемо описати систему n нелінійних, одночасних рівнянь в n невідомих векторів z, векторним рівнянням:

Applying a truncated Taylor's formula in the vicinity of the j-th estimate of the roots, z(j), gives

Застосовуючи скорочену формулу Тейлора в околі J-тих коренів, z(J), дає

To obtain the (j + l)th estimate of the roots, namely z (j+a), we note that these roots should ideally satisfy equation (2.71) precisely:

Для отримання (J + 1) -й оцінки коренів, а саме z (J +a), зауважимо, що ці корені в ідеалі повинні задовольняти рівнянню (2.71), a саме:

Substituting from (2.73) into (2.72) allows us to write down our next estimate for z, z (j+t), as:

Підставляючи з (2.73) в (2.72) дозволяє записати наші наступні оцінки z, z(J + 1), як:

This formula is used to converge iteratively on the true roots of the equation set.

Ця формула використовується для інтерактивного сходження справжніх коренів безлічі рівняння.

To use this result to solve equation (2.70) for the values of the states at the (k + l ) t h timestep, x(k+l), we put

Щоб використовувати цей результат для рішення рівняння (2.70) для значень станів на (к + 1) -го тимчасового кроку, х (к + 1), ми ставимо

Since df(k+t)/dx(k+1) is simply the Jacobian, J, as defined in equation (2.48) at the time instant, t0+(k+l)Δt, we may rewrite (2.76) as:

Так df(k+t)/dx(k+1) є простим Якобі, J, як визначено в рівнянні (2.48) в момент часу, t0+(k+l)Δt, ми можемо переписати (2.76) наступним чином:

Hence we may generate the (j + 1 )th estimate for x(к+1) from the j-th such estimate by:

Отже, ми можемо генерувати (J + 1) ряд на Хк+1 від J-ї такої оцінки по:

A starting value for Xk+l will be guessed (and once the integrations have begun a good value of this vector will be its value at the last timestep, xk) and equation (2.78) can be evaluated repeatedly until some criterion of convergence is satisfied.

Початкове значення для Xk+1 будуть передбачатись (а коли інтегрування починатиметься значення цього вектора буде відповідати його значенню на останньому тимчасовому кроку, Хк) і рівняння (2.78) може виконуватись доти, доки деякі критерієм збіжності не виконаються.

Note that Uk+l and Xk on the right-hand side are invariant throughout all the iterations at each timestep. The Jacobian is evaluated numerically, and strictly this should be done at each iteration. But this is a time-consuming procedure, and in practice the Jacobian is not usually calculated at each iteration, and not even at every timestep. Further time is saved by using sparse matrix techniques to take advantage of the fact that the Jacobian usually possesses many zero elements (cf. equation (2.52) for example). Sparse matrix techniques are similarly used in solving equation (2.78) once the Jacobian has been found. Finally, the integration routine will seek to lengthen the timestep to the maximum extent consistent with a defined accuracy criterion, to take advantage of the strong stability properties of the implicit method.

Зверніть увагу, що Uk+1 і Xk на правій стороні інваріантні у всіх ітерацій на кожному часовому кроці. Якобіан чисельно, і строго, це повинно бути зроблено на кожній ітерації. Але це тривала процедура, і на практиці Якобі як правило, не розраховані на кожній ітерації, і навіть не на кожному часовому кроці. Крім того можна економиться час, скористатися тим, що використовуємо рідкісні матричні методи, а також якобіан зазвичай має багато нульових елементів (див рівняння (2.52), наприклад). Рідкісні матричні прийоми аналогічним використовуються при вирішенні рівняння (2.78) один раз Якобі був знайдений. Нарешті, процедура інтегрування прагнутиме продовжити часовий крок, щоб в максимальному ступені відповідно до певної критерію точності, щоб скористатися сильними властивостями стійкості неявного методу.

As a result of including an implicit, 'stiff' integration algorithm such as the first-order algorithm described in outline above, a simulation package may speed up markedly the execution speed of the simulation as a whole. Several different integration routines, some designed for stiff systems and some not, may well be provided in the simulation package, and it will be for the control engineer to decide which he wishes to use. Often this will be through a process of trial and error, with speed, stability and accuracy as the objectives.

В результаті отримаєм в тому числі неявно, алгоритм «жорсткої» інтеграції, таких як алгоритм першого порядку, описаного в загальних рисах вище, пакет моделювання може помітно прискорити швидкість виконання моделювання в цілому. Кілька різних процедур інтегрування, одні з яких призначені для жорстких систем, і деякі з них не можуть, також бути передбачені в пакеті моделювання, і це буде для інженера-автоматчика, для вибору, який він бажає використовувати. Часто це буде обиратись в процесі проб і помилок, зі швидкістю, стабільності і точності в якості цілей.

2.9 Tackling stiffness in process simulations by modifications to the model
The modeller will not normally wish to tamper with the stiff integration algorithms provided with his modelling package- it would almost always be counterproductive for him to repeat programming carded out by the package designer and already tested to a high degree. Nevertheless, the modeller's physical grasp of the problem can allow him to reduce the stiffness of the equations finally presented for numerical integration. Considering equation (2.44), repeated below

2.9 Вибір жорсткість у процесі моделювання за допомогою модифікацій моделі
Моделювальник як правило, не бажають возитися з жорстким алгоритмом інтегрування за умови, що його пакет для моделювання - це майже завжди контрпродуктивно для його повтору шляхом програмування і вже випробували в високого ступеня жорсткості. Тим не менш, фізичні проблеми при моделюванні можуть дозволити йому зменшити жорсткість рівнянь нарешті представивши його для чисельного інтегрування. Розглядаючи рівняння (2.44), повторюється нижче

it is clear that stiffness may be reduced by increasing the value of Tmin,. This may be done in two ways:
(i) the modeller may decide to increase artificially the minimum time constant so that it is closer to the other time constants dominating what he considers to be the essential behaviour of the model, or
(ii) by assuming that the fastest part of the model responds instantaneously, he may decrease the minimum time constant to zero and thus take that time constant out of consideration.

ясно, що жорсткість може бути зменшена за рахунок збільшення значення Tmin,. Це може бути зроблено двома способами:
(I) моделювальник може прийняти рішення про штучне збільшення мінімальної постійної часу таким чином, що вона ближче до інших домінуючих постійних часу, що буде важливим для поведінка моделі, або
(II), якщо припустити, що найшвидша частина моделі відповідає миттєво, він може зменшити мінімальну постійну часу до нуля, і, таким чином, вважати, що постійна часу знята з розгляду.

It may seem paradoxical but it is nevertheless true that these two opposite courses of action have essentially the same effect on the stiffness ratio, S. In the first case, the effect is to raise the value of the smallest time constant towards the next smallest:

Це може здатися парадоксальним, але це, тим не менш вірно, що ці два протилежні напрямки діяльності мають, по суті той же самий ефект на коефіцієнт жорсткості, S. У першому випадку, ефект полягає в підвищенні значення найменшого постійної часу до наступної меншої:

In the second case, the smallest time constant disappears from consideration, so that the new stiffness ratio is given by the similar formula:

У другому випадку, найменша стала часу зникає з розгляду, так що новий коефіцієнт жорсткості задається аналогічною формулою:

To gain an understanding of the physical significance of these two courses of action, let us refer once more to the tank liquid level system of Section 2.2, working near the operating point set out in Table 2.2. The dominant time constant is that associated with liquid transit time, namely 88.99 s. If our principal concern is merely level control, it will make little difference to the result of the simulation if we increase the values of the time constants associated with valves 1 and 2 to 10 s, say. But by doing so, we reduce the stiffness ratio from 29 to 9, a useful gain. Equally, it will make little difference if we artificially reduce the time constants for valves 1 and 2 to zero. The differential equations of (2.9) and (2.10) are then replaced by the simple algebraics:

Щоб отримати уявлення про фізичну значимісті цих двох напрямків діяльності, звернемося ще раз до системи рівня рідини в резервуарі з розділу 2.2, працюючи поруч з операційною точкою, яка викладена в таблиці 2.2. Домінуючою постійної часу є та, що пов'язана з часом перекачки рідини, а саме 88,99 сек. Якщо наша головна турбота просто регулювання рівня, то мало що зміниться в результаті моделювання, якщо ми збільшуємо значення постійних часу, пов'язаних з клапанами 1 і 2 до 10 с, скажімо. Але тим самим ми зменшуємо коефіцієнт жорсткості з 29 до 9, і отримуєм корисну вигоду. В рівній мірі це мало що змінить, якщо ми штучно знизим постійні часу для клапанів 1 і 2 до нуля. Диференціальні рівняння (2.9) і (2.10) будуть замінені простими алгебраїчними:

The stiffness ratio is reduced to 8, and the simplification has brought the additional bonus of reducing the number of states from 4 to 2.

Співвідношення жорсткості знижується до 8, і спрощення призвело появи додаткового бонусу зменшення кількості станів від 4 до 2.