JP2001290514A - Simulation method of control simulation device - Google Patents

Abstract

PROBLEM TO BE SOLVED: To increase a simulation speed by providing a method for solving a differential equation representing a system fast by an implicit numerical integration method for the computer simulation of a control system. SOLUTION: Such a state variable that the variation quantity of its time differential coefficient is not 0 is specified from the connection relation among the modules constituting the control system and an output equation and only its variation quantity is computed to decrease the computation quantity of elements of a Jacobian matrix, thereby making the simulation fast.

Description

【発明の詳細な説明】DETAILED DESCRIPTION OF THE INVENTION

【０００１】[0001]

【発明の属する技術分野】本発明は、制御系のシミュレ
ーションを行う装置、特に高速なシミュレーションが要
求される制御系のシミュレーション装置におけるシミュ
レーション方法に関する。BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to an apparatus for simulating a control system, and more particularly to a simulation method for a simulation apparatus for a control system which requires a high-speed simulation.

【０００２】[0002]

【従来の技術】自然界の現象の中には、数学的に常微分
方程式で記述できるものが多く、制御の分野において
も、制御対象となるプラントおよびそれを制御するコン
トローラなどを常微分方程式でモデル化（記述）し、そ
の常微分方程式を解析することによりプラントおよびコ
ントローラを含めた系全体の挙動を把握する。そのため
にはモデル化された常微分方程式を解く必要があるが、
しかし解析的に解ける常微分方程式は非常に限られ、現
実的には計算機を用いて数値的に解くという手法（数値
積分法）が採られるのが一般的である。2. Description of the Related Art Many phenomena in the natural world can be mathematically described by ordinary differential equations. In the field of control, a plant to be controlled and a controller for controlling the same are modeled by ordinary differential equations. Then, the behavior of the whole system including the plant and the controller is grasped by analyzing (or describing) the ordinary differential equation. To do this, we need to solve the modeled ODE,
However, the ordinary differential equations that can be solved analytically are extremely limited, and in practice, a method of numerically solving using a computer (numerical integration method) is generally adopted.

【０００３】そのような数値積分法としてこれまで様々
な手法が提案されており、シミュレーションの対象や目
的などによってそれら手法の中から適切な手法を選択し
用いることとなる。主な手法としてはルンゲクッタ法に
代表される陽的積分法、台形法に代表される陰的積分法
などがある。Various methods have been proposed as such a numerical integration method, and an appropriate method is selected from these methods depending on the object or purpose of the simulation. The main methods include an explicit integration method represented by the Runge-Kutta method and an implicit integration method represented by the trapezoidal method.

【０００４】[0004]

【発明が解決しようとする課題】シミュレーションの対
象期間を長く設定する場合、特にシミュレーション対象
が大規模複雑となると、シミュレーションの各時間ステ
ップに要する計算量は膨大なものとなり、シミュレーシ
ョンに要する時間が長くなる。この計算量を減らすため
には時間ステップを大きく設定することが考えられる
が、陽的積分法は、時間ステップを大きく設定すると原
理上計算誤差が大きくなる傾向があり、シミュレーショ
ン結果の信頼性が損なわれる。一方、陰的積分法は、一
般的に時間ステップを大きく設定しても計算誤差が大き
くなりにくいという性質を有しているが、各時間ステッ
プに要する計算量が陽的積分法よりも大きくなり、結果
的にシミュレーションの速度が向上しない場合があると
いう問題を有している。When the simulation target period is set to be long, especially when the simulation target is large-scale and complicated, the amount of calculation required for each time step of the simulation is enormous, and the time required for the simulation is long. Become. To reduce the amount of calculation, it is conceivable to set a large time step.However, in the explicit integration method, if the time step is set large, the calculation error tends to increase in principle, and the reliability of the simulation result is impaired. It is. On the other hand, the implicit integration method generally has the property that the calculation error is unlikely to increase even if a large time step is set, but the amount of calculation required for each time step is larger than the explicit integration method. As a result, there is a problem that the simulation speed may not be improved.

【０００５】[0005]

【課題を解決するための手段】上記課題を解決するため
に成された本発明では、制御系を構成する複数のモジュ
ールをつなぎ合わせることによって１つの制御系を構成
し、各モジュールの動特性を支配する常微分方程式で表
された状態方程式を、非線型連立方程式の求解にニュー
トン法を用いた陰的数値積分法により時間ステップ毎に
解くことによって、毎時間ステップにおける各モジュー
ルの出力値を計算する制御シミュレーション装置のシミ
ュレーション方法において、陰的数値積分法における非
線型連立方程式をニュートン法により求解するに当た
り、各状態変数を微小に変動させ、その変動に対する各
状態変数の時間微係数の変動量を計算し、それらの値か
らニュートン法のヤコビアン行列を算出する際に、各モ
ジュールを、その微小に変動させる状態を有するモジュ
ールの出力を直接入力とするモジュールと、入力を時間
遅れなく出力とするモジュールを１個以上介して入力と
するモジュールと、入力としないモジュールとに、系を
構成するモジュールの接続関係および各モジュールの出
力方程式から選別することで、状態変数の時間微係数の
変動量が０とならない状態変数を特定し、その変動量の
みを計算によって求めることによって、ヤコビアン行列
の各要素の値を計算する際の計算量を減らし、よって数
値積分に要する時間を短縮し、シミュレーションを高速
化する。According to the present invention which has been made to solve the above-mentioned problems, one control system is formed by connecting a plurality of modules constituting a control system, and the dynamic characteristics of each module are determined. Calculate the output value of each module at each time step by solving the equation of state expressed by the governing ordinary differential equation at each time step by solving the nonlinear simultaneous equations using the implicit numerical integration method using Newton's method In solving the nonlinear simultaneous equations in the implicit numerical integration method by the Newton method in the simulation method of the control simulation device, each state variable is minutely changed, and the amount of change of the time derivative of each state variable with respect to the change is calculated. When calculating and calculating the Newton method Jacobian matrix from those values, each module is The modules that constitute a system include a module that directly receives an output of a module that has a state that fluctuates, a module that receives an input via at least one module that outputs an input without time delay, and a module that does not receive an input. By determining from the connection relation of the module and the output equation of each module, the state variable in which the variation of the time derivative of the state variable does not become 0 is specified, and only the variation is obtained by calculation, thereby calculating each element of the Jacobian matrix. Reduces the amount of calculation when calculating the value of, thereby shortening the time required for numerical integration and speeding up the simulation.

【０００６】[0006]

【発明の実施の形態】本実施例の制御シミュレーション
装置を実現する構成を図１に示す。図１で、１０１は中
央演算装置（ＣＰＵ）、１０２はプログラムの格納およ
び演算結果の格納に用いられる常時書き込み読み出し可
能記憶装置（ＲＡＭ）、１０３は計算結果を表示するた
めのディスプレイ装置、１０４はプログラムおよびデー
タを格納するためのハードディスク装置、１０５は１０
１から１０４でを結合するバスである。FIG. 1 shows a configuration for realizing a control simulation apparatus according to this embodiment. In FIG. 1, 101 is a central processing unit (CPU), 102 is a constantly writable and readable storage device (RAM) used for storing programs and calculation results, 103 is a display device for displaying calculation results, and 104 is a display device. Hard disk drive for storing programs and data, 105 is 10
A bus connecting 1 to 104.

【０００７】この構成により、図２のフローチャートに
従って数値積分が行われる。以下図３のようなシミュレ
ーション対象を例にし、陰的積分法として台形法を用い
た場合について説明する。なお、図３の３０１は出力を
ｙ₁ 、状態変数をｘ₁ ，ｘ₂、状態方程式ｄ（ｘ₁ ，ｘ
₂ ）／ｄｔ＝ｆ₁ （ｘ₁ ，ｘ₂ ）、出力方程式をｙ₁＝
ｈ₁ （ｘ₁ ，ｘ₂ ）を動特性として持つモジュール、３
０２は入力をｕ₂ 、状態変数をｘ₃ 、出力をｙ₂ 、状態
方程式をｄｘ₃ ／ｄｔ＝ｆ₂ （ｘ₃ ）＋ｇ₂ （ｕ₂ ）、
出力方程式をｙ₂ ＝ｈ₂ （ｘ₃ ）として持つモジュー
ル、３０３は入力をｕ₃ 、状態変数をｘ₄ ，ｘ₅ 、状態
方程式をｄ（ｘ₄ ，ｘ₅ ）／ｄｔ＝ｆ₃ （ｘ₄ ，ｘ₅ ）
＋ｇ₃ （ｕ₃ ）として持つモジュール、３０４は入力を
ｕ₄ 、出力をｙ₄ 、状態変数をｘ₆ 、状態方程式をｄｘ
₆ ／ｄｔ＝ｆ₄ （ｘ₆ ）＋ｇ₄ （ｕ ₄ ）、出力方程式を
ｙ₄ ＝ｈ₄ （ｘ₆ ）＋ｋ₄ （ｕ₄ ）として持つモジュー
ル、３０５は入力をｕ₅ 、状態変数をｘ₇ 、状態方程式
をｄｘ₇ ／ｄｔ＝ｆ₅ （ｘ₇）＋ｇ₅ （ｕ₅ ）として持
つモジュールである。ここで、３０４は、出力方程式が
その入力ｕ₄ の関数となっているため、入力ｕ₄ が時間
遅れなく出力ｙ₄ に現れるという性質（直達成分）を持
っているモジュールであり、その他のモジュールはその
ような性質は持っていない。このような各モジュールの
直達成分の有無は、直達成分フラグによって各モジュー
ル毎に記憶される。また、モジュール間の接続関係は、
１０２に示すように、各モジュールの次モジュールへの
ポインタによってＲＡＭ上に記憶される。With this configuration, the flowchart of FIG.
Therefore, numerical integration is performed. A simulation as shown in Figure 3 below
Using the trapezoidal method as an implicit integration method
The following describes the case where The output 301 in FIG.
y₁, And state variables x₁, X_Two, The state equation d (x₁, X
_Two) / Dt = f₁(X₁, X_Two), The output equation is y₁=
h₁(X₁, X_Two) As a dynamic characteristic, 3
02 is input u_Two, And state variables x_Three, Output y_Two,Status
Dx equation_Three/ Dt = f_Two(X_Three) + G_Two(U_Two),
Output equation is y_Two= H_Two(X_Three) Module
303, input u_Three, And state variables x_Four, X_Five,Status
The equation is d (x_Four, X_Five) / Dt = f_Three(X_Four, X_Five)
+ G_Three(U_Three)), 304 is an input
u_Four, Output y_Four, And state variables x₆And the state equation is dx
₆/ Dt = f_Four(X₆) + G_Four(U _Four), The output equation
y_Four= H_Four(X₆) + K_Four(U_Four) Module
305 input u_Five, And state variables x₇, Equation of state
To dx₇/ Dt = f_Five(X₇) + G_Five(U_Five)
One module. Where 304 is the output equation
That input u_FourInput u_FourIs time
Output y without delay_FourHas the property of appearing in
The other modules are
It does not have such properties. Of each such module
The direct achievement flag is set for each module by the direct achievement flag.
Is stored for each file. The connection relationship between modules is
As shown in 102, each module is assigned to the next module.
It is stored on the RAM by the pointer.

【０００８】またシミュレーションの現在時刻をｔ、時
間ステップ幅をｈ、状態変数をまとめたベクトルをＸ
（＝（ｘ₁ ，ｘ₂ ，ｘ₃ ，ｘ₄ ，ｘ₅ ，ｘ₆ ，
ｘ₇ ））、関数Ｆ（Ｘ）をＦ（Ｘ）＝（ｆ₁ （ｘ₁ ，ｘ
₂ ）ｆ₂ （ｘ₃ ）ｆ₃ （ｘ₄ ，ｘ₅ ）ｆ ₄ （ｘ₆ ）ｆ₅
（ｘ₇ ））、時刻ｔにおけるＸをＸ_t 、時刻ｔ＋ｈにお
けるＸをＸ_t+h と表す。[0008] Also, the current time of the simulation is t, hour
H is the interval step width, and X is a vector summarizing the state variables.
(= (X₁, X_Two, X_Three, X_Four, X_Five, X₆,
x₇)) And the function F (X) is expressed as F (X) = (f₁(X₁, X
_Two) F_Two(X_Three) F_Three(X_Four, X_Five) F _Four(X₆) F_Five
(X₇)), X at time t is X_tAt time t + h
X to X_{t + h}It expresses.

【０００９】まず、Ｓ２０１において初期状態Ｘ₀ を０
に設定し、Ｓ２０２において現在時刻ｔを０に設定す
る。また時間ステップｈは操作者が任意に与えるものと
する。これらは１０２に示すようにＲＡＭに記憶され
る。Ｓ２０３において時刻ｔから時間ステップｈ後の状
態Ｘ_t+h を非線型連立方程式 （Ｘ_t ＋Ｘ_t+h ）／ｈ＝（Ｆ（Ｘ_t ）＋Ｆ（Ｘ_t+h ））／２ （１） を満たすようなＸ_t+h を図４に示すステップに従ってニ
ュートン法を用いて求める。そしてＳ２０４において現
在時刻をｈだけ進め、Ｓ２０５において終了判断が行わ
れ、まだ終了していない場合はＳ２０３に戻る。First, in step S201, the initial state X _{0 is set} to 0.
, And the current time t is set to 0 in S202. The time step h is arbitrarily given by the operator. These are stored in RAM as shown at 102. In S203, the state _{Xt +} h after the time step h from the time t is converted into the nonlinear simultaneous equation ( _Xt + _{Xt + h} ) / h = (F ( _Xt ) + F ( _{Xt + h} )) / 2 (1) _{Xt + h} that satisfies is satisfied by using the Newton method according to the steps shown in FIG. Then, in S204, the current time is advanced by h. In S205, an end determination is made, and if it has not ended, the process returns to S203.

【００１０】次に図４のニュートン法について述べる。
式（１）を満たすＸ_t+h は図４に示すステップに従って
ニュートン法を用いて求める。Ｓ４０１では、式（１）
の解Ｘ_t+h の暫定解Ｘ_t+h ⁽⁰⁾ を適当に定め、Ｓ４０２
では、そのｉ行ｊ列（ｉ，ｊ＝１７）の要素が Ｊ（ｉ，ｊ）＝δ_ij＋０．５×ｈ×（Δ（ｄｘ_i ／ｄｔ）／Δｘ_j ） （２） で定義されるヤコビアン行列（以下Ｊ行列と呼ぶ）を計
算する。ここで、式（２）におけるδ_ijはクロネッカー
のデルタ δ_ij＝１（ｉ＝ｊの場合） ０（ｉ≠ｊの場合） （３） を表し、Δ（ｄｘ_i ／ｄｔ）／Δｘ_j は、状態変数ｘ_j
がΔｘ_j だけ微小変動したときの状態変数ｘ_i の時間微
係数の微小変動量を表しており、この値は各モジュール
の状態方程式および出力方程式から計算される。例え
ば、ｊ＝１，ｉ＝３の場合は、ｘ₁ がΔｘ₁ だけ微小変
動したときのモジュール３０１の出力ｙ₁（＝ｕ₂ ）を
計算し、ｕ₂ からモジュール３０２の状態変数ｘ₃ の時
間微係数の値を求め、予め計算しておいたｘ₁ が微小変
動する前のｘ₃ の時間微係数の値からの差分を計算する
ことで求められる。次に、Ｓ４０３でＪ行列を用いて解
の更新方向ｄを線形連立方程式 Ｊｄ＝−Ｆ（Ｘ_t+h ） （４） を解き求める。次に、Ｓ４０４において Ｘ_t+h ⁽ⁿ⁺¹⁾ ←Ｘ_t+h ⁽ⁿ⁾ ＋μｄ （５） と更新する。そしてＳ４０５において解の収束を判定
し、収束していない場合Ｓ４０２に戻り、収束している
場合は終了する。Next, the Newton method of FIG. 4 will be described.
X _{t + h} that satisfies the equation (1) is obtained by using the Newton method according to the steps shown in FIG. In S401, equation (1)
Define the solution X _{t + h} of interim solution X _{t + h} ⁽⁰⁾ appropriately, S402
So the elements of the row i and column j (i, j = 17) is defined by the J (i, j) = δ ij + 0.5 × h × (Δ (dx i / dt) / Δx j) (2) A Jacobian matrix (hereinafter referred to as a J matrix) is calculated. Here, [delta] _ij in equation (2) (in the case of i = j) delta [delta] _ij = 1 Kronecker (for i ≠ j) 0 represents a _{(3), Δ (dx i} / dt) / Δx j is , State variable x _j
Represents the minute variation of the time derivative of the state variable x _i when the value fluctuates by Δx _j , and this value is calculated from the state equation and output equation of each module. For example, when j = 1 and i = 3, the output y ₁ (= u ₂ ) of the module 301 when x ₁ slightly fluctuates by Δx ₁ is calculated, and the state variable x ₃ of the module 302 is calculated from u ₂ . obtains the value of the time derivative, obtained by x ₁ which had been previously calculated to calculate the difference from the value of the time derivative of the previous x ₃ to slight change. Next, in S403, the solution update direction d is solved by using the J matrix to solve the linear simultaneous equation Jd = −F (X _{t + h} ) (4). Next, in S404, Xt _{+ h} ^{(n + 1)} ← Xt _{+ h} ⁽ⁿ⁾ + μd (5) is updated. Then, in step S405, the convergence of the solution is determined. If the convergence has not occurred, the process returns to step S402, and if the convergence has occurred, the process ends.

【００１１】Ｓ４０２では、Ｊ行列を計算する必要があ
るが、従来技術では全てのｉ，ｊについて式（２）を計
算しＪ行列を計算していたが、シミュレーション対象が
大規模複雑になるに従い状態数が膨大になり、式（２）
の計算量が膨大になり、結果シミュレーションの速度が
著しく低くなる。例えば図３では、状態変数が７個ある
ので式（２）の計算を７² ＝４９回行う必要がある。In S402, it is necessary to calculate the J matrix. In the prior art, the formula (2) is calculated for all i and j to calculate the J matrix. The number of states becomes enormous, and equation (2)
The amount of calculation becomes huge, and as a result, the speed of the simulation becomes extremely low. For example, in FIG. 3, since there are seven state variables, the calculation of equation (2) needs to be performed 7 ² = 49 times.

【００１２】そこで本発明では、Ｊ行列を計算する際に
全てのｉ，ｊについて式（２）を計算するのではなく、
Ｊ行列のうち０となる要素を前もって特定しておき、そ
の要素に関しては式（２）の計算を省く事でＪ行列を計
算する際の全体の計算量を低減させる。一般的に、シミ
ュレーション装置におけるＪ行列は０要素が多いため、
大幅な計算量の低減が期待される。Therefore, in the present invention, when calculating the J matrix, instead of calculating equation (2) for all i and j,
An element that becomes 0 in the J matrix is specified in advance, and the calculation amount of the J matrix is reduced by omitting the calculation of Expression (2) for the element, thereby reducing the total amount of calculation when calculating the J matrix. In general, the J matrix in a simulation device has many 0 elements,
A large reduction in the amount of calculation is expected.

【００１３】本発明でのＪ行列の計算は図５に示すステ
ップに従って計算する。Ｓ５０１において、変動させる
状態変数ｘ_j （ｊ＝１〜７）を任意に選択する。そして
Ｓ５０２においてそのｘ_j の変動幅Δｘ_j をパラメータ
により適当に設定する。Ｓ５０３においてｘ_j の変動の
影響を求める状態変数ｘ_i を任意に選択する。The calculation of the J matrix in the present invention is performed according to the steps shown in FIG. In S501, a state variable x _j (j = 1 to 7) to be changed is arbitrarily selected. And appropriately setting the variation range [Delta] x _j of the x _j by the parameter in S502. Arbitrarily select a state variable x _i to determine the effects of changes in x _j in S503.

【００１４】Ｓ５０４においてｘ_j を有するモジュール
の出力がｘ_i を有するモジュールの入力に直接接続され
ているかを、１０２のＲＡＭ上に記憶した各モジュール
の次モジュールポインタを順次辿って判断し、直接接続
されていない場合Ｓ５０５で途中１つ以上のモジュール
を経由して間接的に接続しているかを、同じく１０２の
ＲＡＭ上に記憶された各モジュールの次モジュールへの
ポインタを辿ることで判断する。直接接続されている場
合、および間接的に接続されかつ途中経由しているモジ
ュールがすべて直達成分（入力信号が微分要素もしくは
積分要素を経由せず出力信号に現れる成分）を有してい
る場合、Ｓ５０７において式（２）に従ってＪ行列のｉ
行ｊ列の成分を計算する。In step S504, it is determined whether the output of the module having x _j is directly connected to the input of the module having x _i by sequentially following the next module pointer of each module stored in the RAM of 102, and is directly connected. If not, it is determined in step S505 whether the connection is indirectly via one or more modules in the middle by tracing the pointer to the next module of each module stored in the RAM 102. When directly connected, and when all the modules connected indirectly and passing therethrough have direct achievements (components in which the input signal does not pass through the differential element or the integral element and appears in the output signal), In S507, i of the J matrix according to the equation (2)
Calculate the components in row j column.

【００１５】ここで、直達成分の有無は１０２のＲＡＭ
上に記憶された各モジュールの直達成分存在フラグによ
って判断する。それ以外の場合は式（２）の計算を行う
までもなくＪ行列のｉ行ｊ列の成分は０となるため、式
（２）の計算を省略する事ができ、計算量を削減する事
ができる。また、ｘ_i とｘ_j が同じモジュールに含まれ
ている場合もＪ行列のｉ行ｊ列の成分は０となり、ｉ＝
ｊの場合は式（２）より１となる。例えば、ｊ＝１の場
合、状態変数ｘ_j を有しているモジュールは３０１であ
り、３０１の出力を直接入力としているのは３０２と３
０４であり、また３０５は直達項をもつモジュール３０
４を経由して３０１に接続しているので、状態変数ｘ_j
が微少変動する際に計算すべき時間微係数は、Δ（ｄｘ
₃ ／ｄｔ）／Δｘ₁ 、Δ（ｄｘ₆ ／ｄｔ）／Δｘ₁ 、Δ
（ｄｘ₇ ／ｄｔ）／Δｘ₁ であり、これらの計算結果と
式（２）からＪ（３，１），Ｊ（６，１），Ｊ（７，
１）を計算する。また、状態変数ｘ₂ はｘ₁ と同じモジ
ュールに含まれているのでＪ（２，１）＝０となり、ま
たＪ（１，１）は式（２）から１となる。状態変数
ｘ ₄ ，ｘ₅ を有するモジュールは、３０１の出力を直接
入力としておらず、また直達成分をもつモジュールを経
由しても接続されていないためＪ（４，１），Ｊ（５，
１）は０となる。次に、Ｓ５０９において全てのｘ_j お
よび全てのｘ_j について行ったか判断を行い、行った場
合終了し、それ以外の場合次のｘ_i もしくはｘ_j を選択
し、同様の処理を全てのｘ_i およびｘ_j について繰り返
す。Here, the presence or absence of the direct achievement is determined by 102 RAMs.
According to the direct achievement presence flag of each module stored above
I judge. Otherwise, calculate equation (2)
Since the component at the i-th row and the j-th column of the J matrix becomes 0 soon, the expression
The calculation of (2) can be omitted, and the calculation amount can be reduced.
Can be. Also, x_iAnd x_jAre included in the same module
Also, the component of the i-th row and the j-th column of the J matrix becomes 0, and i =
In the case of j, it becomes 1 according to equation (2). For example, if j = 1
The state variable x_jThe module having
The reason why the output of 301 is directly input is 302 and 3
04 and a module 305 having a direct term
4, the state variable x_j
Is small when Δ fluctuates slightly, Δ (dx
_Three/ Dt) / Δx₁, Δ (dx₆/ Dt) / Δx₁, Δ
(Dx₇/ Dt) / Δx₁And these calculation results and
From equation (2), J (3,1), J (6,1), J (7,
Calculate 1). Also, the state variable x_TwoIs x₁Same moji as
J (2,1) = 0 because it is included in the
J (1,1) becomes 1 from Expression (2). State variables
x _Four, X_FiveModule directly outputs the output of 301
The module that is not input and has direct achievement
J (4,1), J (5,
1) becomes 0. Next, in step S509, all x_jYou
And all x_jJudge whether you went about
Ends, otherwise the next x_iOr x_jchoose
And the same processing is performed for all x_iAnd x_jRepeat about
You.

【００１６】以上の方法により、図３で示されるシミュ
レーション対象では、Ｊ（１，１），Ｊ（２，２），Ｊ
（３，３），Ｊ（４，４），Ｊ（５，５），Ｊ（６，
６），Ｊ（７，７）が１、Ｊ（３，１），Ｊ（６，
１），Ｊ（７，１），Ｊ（３，２），Ｊ（６，２），Ｊ
（７，２），Ｊ（４，３），Ｊ（５，３），Ｊ（７，
６）について式（２）を用いて値を計算、残りの残り３
２要素は式（２）を計算するまでもなく０と特定でき、
式（２）の計算回数を１６回に低減する事ができる。According to the above method, in the simulation object shown in FIG. 3, J (1,1), J (2,2), J (1,2)
(3,3), J (4,4), J (5,5), J (6,
6), J (7,7) is 1, J (3,1), J (6,
1), J (7, 1), J (3, 2), J (6, 2), J
(7,2), J (4,3), J (5,3), J (7,
For 6), calculate the value using equation (2), and calculate the remaining 3
Two elements can be specified as 0 without calculating equation (2),
The number of calculations of the equation (2) can be reduced to 16 times.

【００１７】[0017]

【発明の効果】本発明により、シミュレーション対象で
あるモデルの特性を用いて、陰的数値積分法を用いた制
御シミュレーションの実行に要する計算量を低減するこ
とができ、これまで知られている制御シミュレーション
装置よりも高速な制御シミュレーションが可能となるこ
とから、特に大規模複雑な制御シミュレーションを高速
に実行する用途に有用される。According to the present invention, the amount of calculation required to execute a control simulation using an implicit numerical integration method can be reduced using characteristics of a model to be simulated. Since the control simulation can be performed at a higher speed than that of the simulation device, it is particularly useful for the purpose of executing a large-scale and complicated control simulation at a high speed.

【図面の簡単な説明】[Brief description of the drawings]

【図１】本実施例の装置構成を示す図である。FIG. 1 is a diagram showing a device configuration of the present embodiment.

【図２】本実施例の制御シミュレーションの手順を示す
図である。FIG. 2 is a diagram illustrating a procedure of a control simulation according to the present embodiment.

【図３】本実施例の制御シミュレーションのモデル例を
表す図である。FIG. 3 is a diagram illustrating a model example of a control simulation according to the present embodiment.

【図４】本実施例の制御シミュレーションにおけるニュ
ートン法を用いた非線型連立方程式の求解の手順を示す
図である。FIG. 4 is a diagram showing a procedure for solving a non-linear system of equations using the Newton method in the control simulation of the present embodiment.

【図５】本実施例におけるニュートン法のヤコビアン行
列を計算する手順を示す図である。FIG. 5 is a diagram showing a procedure for calculating a Jacobian matrix of the Newton method in the embodiment.

【符号の説明】[Explanation of symbols]

１０１…ＣＰＵ １０２…ＲＡＭ １０３…ＣＲＴ １０４…ＨＤＤ １０５…バス Ｓ２０１〜Ｓ２０５…本実施例の制御シミュレーション
の手順である。 Ｓ３０１〜Ｓ３０５…本実施例の制御シミュレーション
のモデル例における各モジュールである。 Ｓ４０１〜Ｓ４０５…本実施例の制御シミュレーション
におけるニュートン法の手順である。 Ｓ５０１〜Ｓ５０９…本実施例におけるニュートン法の
ヤコビアン行列を計算する手順である。101: CPU 102: RAM 103: CRT 104: HDD 105: Buses S201 to S205: Control simulation procedures of the present embodiment. S301 to S305 are modules in the control simulation model example of the present embodiment. S401 to S405 are the procedures of the Newton method in the control simulation of the present embodiment. S501 to S509 are the procedures for calculating the Jacobian matrix of the Newton method in this embodiment.

Claims (1)

【特許請求の範囲】[Claims] 【請求項１】 制御系を構成する複数のモジュールをつ
なぎ合わせることによって１つの制御系を構成し、各モ
ジュールの動特性を支配する常微分方程式で表された状
態方程式を、非線型連立方程式の求解にニュートン法を
用いた陰的数値積分法により時間ステップ毎に解くこと
によって、毎時間ステップにおける各モジュールの出力
値を計算する制御シミュレーション装置におけるシミュ
レーション方法において、前記ニュートン法により求解
するに当たり、前記状態変数を微小に変動させ、その変
動に対する各状態変数の時間微係数の変動量を計算し、
それらの値からニュートン法のヤコビアン行列を算出す
る際に、各モジュールから、微小に変動させる状態変数
を有するモジュールの出力を直接入力とするモジュール
と、入力を時間遅れなく出力とするモジュールを１つ以
上介して入力とするモジュールとを微少に変動させる状
態変数毎に選別して、前記選別したモジュールに含まれ
る状態変数だけの時間微係数の変動量を計算することに
より、ヤコビアン行列の各要素の値を計算する際の計算
量を減らし、シミュレーションを高速化することを特徴
とする制御シミュレーション装置におけるシミュレーシ
ョン方法。1. A control system is constructed by connecting a plurality of modules constituting a control system, and a state equation expressed by an ordinary differential equation governing dynamic characteristics of each module is converted into a nonlinear simultaneous equation. By solving for each time step by an implicit numerical integration method using Newton's method for solving, in a simulation method in a control simulation device for calculating an output value of each module in each time step, when solving by the Newton's method, The state variable is fluctuated minutely, and the fluctuation amount of the time derivative of each state variable with respect to the fluctuation is calculated,
When calculating the Jacobian matrix of the Newton's method from these values, one module directly receives the output of a module having a state variable that fluctuates minutely, and one module outputs the input without time delay from each module. The module to be input through the above is selected for each state variable that fluctuates minutely, and the amount of change of the time derivative of only the state variable included in the selected module is calculated, whereby each element of the Jacobian matrix is calculated. A simulation method in a control simulation device, characterized in that a calculation amount for calculating a value is reduced and a simulation is speeded up.