このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。
多目的最適化を用いた核燃料処分計画
この例では、いくつかの整数制約を持つ大規模な非線形多目的問題を定式化して解決する方法を示します。この問題はMontonen、Ranta、Mäkelä [1]から改変したものです。目標は、使用済み核燃料を処分することであり、コストを最小限に抑え、使用済み核燃料集合体を原子炉から取り出してから埋め立てるまでの時間を最小限に抑え、一度に保管される使用済み核燃料集合体の数を最小限に抑えることを目指しています。この問題は複数期間の計画問題であり、各期間は 5 年間です。
モデルの概要
原子炉は廃棄物を生み出し、長期にわたって処分するためには埋め立てなければなりません。これらの廃棄物は、使用済み核燃料を含む燃料棒集合体です。原子炉から取り出されたときは、アセンブリは高温で放射能を帯びていますが、徐々にその放射能は低下します。周期 t では、放射線レベルはパラメータ ベクトル u を持つ二重指数関数的減衰関数によって与えられます。
.
各アセンブリは、中間貯蔵施設に移送できるほど十分に冷えるまで原子炉建屋内のプールに保管され、そこで再びプールに保管されます。アセンブリがさらに冷却されると、他のアセンブリとともに銅鉄製の容器に封入され、廃棄トンネルに埋められます。すべての廃棄トンネルは中央トンネルに接続されています。
この図は、原子炉から最終処分までの核燃料集合体の段階を示しています。
問題の変数は、各時間単位が 5 年を表すスケジュールに関連しています。期間は 1 から始まります。
モデルに関連付けられた定数
Z は、燃料集合体が原子炉から取り出される最後の期間です。削除期間は 1:Z です。[1]ではZ = 11である。各期間は 5 年なので、燃料集合体が取り外される最後の期間は 55 時点になります。
N は、キャニスターが埋められる最後の期間です。埋葬期間は1:Nです。[1]ではN = 19である。各期間は 5 年なので、キャニスターを埋めることができる最後の期間は 95 の時点です。
Z = 11; N = 19;
a は、最初の廃棄前の最後の除去期間です。
b は、最後の除去が行われる廃棄期間です。
a = 5; b = 6;
R は、アセンブリを格納するための最小期間数です。
R = 4;
K は、1 つのキャニスターに収まるアセンブリの最大数です。
K = 4;
T は、1 期間に廃棄されるキャニスターの最小数です。
T = 50;
U は、1 期間に廃棄されるキャニスターの最大数です。
U = 500;
Q は廃棄トンネルの長さ（メートル単位）です。
Q = 350;
M(i) は、時刻 i に削除されたアセンブリの数です。
M = 300 - 60*(-1).^(1:Z); % 360 for odd indices, 240 for evenA(i,j) は、期間 j における除去 i からのアセンブリの保管時間です (i <= Z および j <= N )。
A = zeros(Z,N); for i = 1:Z for j = i:N A(i, j) = j - i; end end
p(i,j) は、期間 j における除去 i からのアセンブリの崩壊熱出力 (ワット単位) です (i <= Z および j <= N )。
% The following parameters fit P_{i,j} of Table A2 from [1] to within 1 in each % entry (fractional error <= 1/250). u = [503 0.1346 260 0.0231]; myfun = @(d)round(u(1)*exp(-u(2)*d) + u(3)*exp(-u(4)*d)); PP = myfun(1:N); pij = zeros(Z,N); for i = 2:Z for j = 1:i pij(i,j) = 1e3; % Dummy values because j >= i does not occur. end end for i=1:Z pij(i,i:end) = PP(1:(N-i+1)); % Same decay profile for all removal times end
pmaxup はキャニスターの平均出力の上限であり、pmaxlow は下限です。
pmaxup = 1830; pmaxlow = 1300;
dDTup は処分トンネル間の距離の上限であり、dDTlow は下限です。
dDTup = 50; dDTlow = 25;
dCAup は廃棄トンネル内のキャニスター間の距離の上限であり、dCAlow は下限です。
dCAup = 15; dCAlow = 6;
[1]では手術に伴う費用は示されていない。この例では、次の値が想定されています。
Casは、期間ごとの 1 つのアセンブリの保管コストです。Cisは、中間保管における期間ごとの 1 つのアセンブリの保管コストです。Cspは、アセンブリあたりの保管場所のコストです。Ccaはキャニスターのコストです。Cefは、期間あたりのカプセル化施設の運用コストです。Cdtは廃棄トンネルの 1 メートルあたりのコストです。Cctは中央トンネルの 1 メートルあたりのコストです。
Cas = 50; Cis = 60; Csp = 10; Cca = 1200; Cef = 300; Cdt = 3000; Cct = 5000;
問題の最適化変数
MATLAB® の問題を作成するには、問題ベースのアプローチを使用します。ほとんどの量に対して連続変数を定義します。
この図は、アセンブリの動きに関連する変数を示しています。
x(i,j) は、時刻 i に削除され、時刻 j、i <= Z、および j <= N に処分されたアセンブリの数です。
x = optimvar("x",Z,N,LowerBound=0,UpperBound=U*K);z(i,j) は、 i の時点で取り外され、 j、 i <= Z、 j <= N の時点で保管されていたアセンブリの数です。
z = optimvar("z",Z,N,LowerBound=0,UpperBound=U*K*N/2);y(j) は、時刻 j <= N に廃棄されたキャニスターの数です。
y = optimvar("y",N,LowerBound=0,UpperBound=U);pmax はキャニスターの最大平均出力です。
pmax = optimvar("pmax",LowerBound=pmaxlow,UpperBound=pmaxup);この図は、廃棄トンネルに関連する数量を示しています。
dDT は隣接する廃棄トンネル間の距離です。
dDT = optimvar("dDT",LowerBound=dDTlow,UpperBound=dDTup);dCA は廃棄トンネル内の隣接するキャニスター間の距離です。この距離は、この例の「問題の制約」セクションで、関数 g を使用して後で計算します。
この数字はカプセル化時間に関連しています。
次の変数を、下限が 0、上限が 1 の整数型バイナリ変数として指定します。
ej(j) は、期間 j <= N 中にカプセル化機能が動作していることを示します。
ej = optimvar("ej",N,Type="integer",LowerBound=0,UpperBound=1);
ejON(j) は、期間 j <= N の初めにカプセル化が開始されることを示します。
ejON = optimvar("ejON",N,Type="integer",LowerBound=0,UpperBound=1);
ejOFF(j) は、期間 j <= N の初めにカプセル化が終了することを示します。
ejOFF = optimvar("ejOFF",N,Type="integer",LowerBound=0,UpperBound=1);
この数字は、アセンブリを破棄できる時間、rij = 1 に関係します。これらの時間は ejON = 1 に始まり、 sij = 1 に終わります。
sij(i,j) は、時間 i に削除されたアセンブリが、時間 j、i <= Z、および j <= N の開始時から廃棄できなくなることを示します。
sij = optimvar("sij",Z,N,Type="integer",LowerBound=0,UpperBound=1);
rij(i,j) は、時刻 i に削除されたアセンブリが時刻 j、i <= Z、および j <= N に処分できることを示します。
rij = optimvar("rij",Z,N,Type="integer",LowerBound=0,UpperBound=1);
すべての最適化変数と問題パラメータが定義されました。
問題の制約
目的と制約を保持する最適化問題を作成します。
prob = optimproblem;
制約数は[1]の式と一致する。最初の 3 つの制約は、一時保管中のアセンブリの数に関係します。
jnot1 = 2:N; prob.Constraints.cons10 = z(:,1) - M(:) + x(:,1) == 0; prob.Constraints.cons11 = z(:,jnot1) - z(:,(jnot1 - 1)) + x(:,jnot1) == 0; prob.Constraints.cons12 = z(:,N) == 0;
すべてのアセンブリが一度破棄されるという制約を設定します。
prob.Constraints.cons13 = sum(sij,2) == 1;
次の制約を設定して変数 rij を定義します。
cons15 = optimconstr(Z,N); cons15(:,1) = rij(:,1) == -sij(:,1) + ejON(1); % equation 14 cons15(:,jnot1) = rij(:,jnot1) == ... rij(:,jnot1-1) - sij(:,jnot1) + repmat(ejON(jnot1)',Z,1); % equation 15 prob.Constraints.cons15 = cons15;
カプセル化機能が稼働しているときのみ廃棄が行われるという制約を設定します。
cons16 = rij <= repmat(ej', Z, 1); prob.Constraints.cons16 = cons16;
生産能力を超えないという制約を指定します。
prob.Constraints.cons17 = x <= U*K*rij;
次の制約は、廃棄する前にアセンブリが十分に冷却されることを強制します。
prob.Constraints.cons18 = x.*(A - R) >= 0;
次の制約はカプセル化機能に関連しています。これらの制約により、施設のオン/オフは1回のみとなり、すべてのキャニスターが1回の実行でカプセル化されます。
prob.Constraints.cons19 = sum(ejON) == 1; prob.Constraints.cons20 = sum(ejOFF) == 1;
次の制約を設定して変数 ej を定義します。
cons21 = optimconstr(N); cons21(1) = ej(1) == ejON(1) - ejOFF(1); % equation 21 cons21(jnot1) = ej(jnot1) == ... ej(jnot1 - 1) + ejON(jnot1) - ejOFF(jnot1); % equation 22 prob.Constraints.cons21 = cons21;
キャニスターの数が廃棄に十分であり、生産能力を超えず、最小生産制約に従うという制約を設定します。
prob.Constraints.cons23 = y' >= (1/K)*sum(x,1); prob.Constraints.cons24 = y <= U*ej; jnotN = 1:(N-1); prob.Constraints.cons25 = y(jnotN) >= T*(ej(jnotN) - ejOFF(jnotN + 1));
処分施設については、キャニスターの熱出力が許容できるという制約を設定する。
prob.Constraints.cons26 = sum(pij.*x,1) <= pmax*y';
埋設されたキャニスター間の距離に非線形制約を指定します。この関数は区分線形であり、最適化式ではサポートされていない演算である max 関数を使用して定義されます。したがって、制約を prob に配置するには、 fcn2optimexpr を使用します。
g = @(pmax,dDT)max([-2.26911*dDT + 0.00675*pmax + 54.5288,... -0.05833*dDT + 0.00596*pmax - 0.727083,... -0.14*dDT + 0.17701*pmax - 350.651]); dCA = fcn2optimexpr(@(pmax,dDT)g(pmax,dDT),pmax,dDT); prob.Constraints.cons29a = dCA >= dCAlow; prob.Constraints.cons29b = dCA <= dCAup;
コスト目標
この多目的問題の最初の目的はコストであり、コストには 7 つの要素があります。
cost = optimexpr(7,1);
1.アセンブリの保管コスト。このコストは、単位時間あたりのコストと各アセンブリが保管される時間の長さを乗じた合計です。
cost(1) = Cas*sum(A.*x,"all");2.一時保管にかかる費用。このコストは、ejOFF の 1 つのコンポーネントである 1 に対して j*ejOFF(j) です。
cost(2) = Cis*max(ejOFF(1)–1,2*ejOFF(2)–1,3*ejOFF(3)–1,...,N*ejOFF(N)–1).
この式を簡単に表すと、新しい最適化変数uのcost(2) = Cis*uと制約
ucons = u >= ((1:N)'.*ejOFF) – 1.
u = optimvar("u",LowerBound=0);
cost(2) = Cis*u;
ucons = u >= ((1:N)'.*ejOFF) - 1;
prob.Constraints.ucons = ucons;3.アセンブリストレージの位置のコスト。cost(3)はCsp*v1で表すことができます。v1は新しい最適化変数であり、制約とともに
v1 >= sum(M)
1 <= j <= Nごとにv1 >= sum_{i=1}^j z(i,j)
b <= j <= N ごとに v1 >= sum_{i=1}^Z z(i,j) 。
これらのコストと関連する制約を作成します。
v1 = optimvar("v1",LowerBound=0); cost(3) = Csp*v1; % Include the three v1 constraints given below. v1consa = v1 >= sum(M); bmin1 = 1:(b-1); v1consb = optimconstr(b-1); for j=bmin1 v1consb(j) = sum(z(1:a+j,j)) <= v1; end ell = b:N; v1consc = sum(z(:,ell),1) <= v1; prob.Constraints.v1consa = v1consa; prob.Constraints.v1consb = v1consb; prob.Constraints.v1consc = v1consc;
4.キャニスターのコスト。このコストは、キャニスター 1 個あたりのコストと、埋められたキャニスターの総数を掛けたものです。
cost(4) = Cca*sum(y);
5.カプセル化施設の運営コスト。このコストは、単位時間あたりのコストに施設の稼働時間を掛けたものです。
cost(5) = Cef*sum(ej);
6.廃棄トンネルの費用。これは、単位長さあたりのコストにキャニスター間の長さを掛け、さらに埋められたキャニスターの総数を掛けたものです。
cost(6) = Cdt*dCA*sum(y);
7.中央トンネルの費用。このコストは、単位長さあたりのコストに中央トンネルの必要な長さを掛けたものです。1 つの処分トンネルに埋められるキャニスターの数は、その長さ Q をキャニスター間の距離 dCA で割った数です。中央トンネルの長さは、埋設されたキャニスターの数sum(y)に比例し、Q/dCAに反比例し、コストはCctに比例します。
cost(7) = Cct/Q*dDT*dCA*sum(y);
総コストは 7 つのコスト要素の合計です。総コストのスケールを他の目的のスケールと一致するように変更するには、合計の対数を取ります。
prob.Objective.cost = log(sum(cost));
安全目標
この問題には、安全性に関連する 2 つの目的があります。目標 2 (safety1 と呼ばれます) は、すべての削除における最大保存時間を最小限に抑えることを試みます。目標 3 は safety2 と呼ばれ、できるだけ早く廃棄を停止しようとします。このスクリプトの最後に表示されるヘルパー関数 max1 と max2 を使用して、これら 2 つの目的を定義します。
prob.Objective.safety1 = fcn2optimexpr(@max1,A,sij); % Minimize maximum storage time, objective (2) in [1] prob.Objective.safety2 = fcn2optimexpr(@max2,ejOFF); % Stop disposal as early as possible, objective (5) in [1]
オプションの設定
ソルバーが進行するにつれて、パレート図のオプションを設定します。この問題には 900 を超える変数があるため、デフォルトよりも大きい 500 の母集団サイズを使用するようにオプションを設定します。また、問題にはバイナリ変数が含まれているため、mutationgaussian 変異関数を使用します。この変異関数は、バイナリ変数のデフォルトの mutationpower よりも効果的に機能します。
opts = optimoptions("gamultiobj",PlotFcn="gaplotpareto",PopulationSize=5e2,... MutationFcn=@mutationgaussian);
問題の実行
問題の定式化が完了し、この多目的問題のオプションが設定されました。問題を実行します。
rng default % For reproducibility [sol,fval,exitflag,output] = solve(prob,Options=opts);
Solving problem using gamultiobj. gamultiobj stopped because the average change in the spread of Pareto solutions is less than options.FunctionTolerance.
xlabel("Cost") ylabel("Safety 1") zlabel("Safety 2")
パレート図は、コストと安全性 1 の間に明確なトレードオフがあることを示しています。実際のコストは表示されている金額の指数関数であるため、トレードオフは図示されているよりも厳しくなります。
解の検証
gamultiobj は、さまざまな適合関数値を持ついくつかの実行可能なソリューションを見つけます。ソリューションに関連付けられた制御変数を見つけるには、データ ヒントを使用します。
データヒントを有効にした後、左上のソリューションをクリックします。
選択されたポイントのインデックスは 4 です。この点に関連付けられた変数は sol(4) にあります。
このソリューションに関連付けられている x 変数を調べます。x(i,j) は、時間 i に削除され、時間 j に処分されることを思い出してください。
disp(sol(4).x)
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 360.0000 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 240.0000 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 265.4535 94.5465 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 240.0000 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 360.0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 240.0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 254.5465 105.4535 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 240.0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 360.0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 240.0000 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 360.0000 0
明らかに、x スケジュールは整数値に制限されません。除去量M(:)と比較して、処分に関するxスケジュールの合計を表示します。
disp([sum(sol(4).x,2),M(:)])
360 360 240 240 360 360 240 240 360 360 240 240 360 360 240 240 360 360 240 240 360 360
x スケジュールはすべての削除を考慮します。
カプセル化施設の稼働時間は何時ですか？
disp(sol(4).ej')
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0
カプセル化は 17 回目と 18 回目に実行されます。
廃棄トンネル間の距離はどれくらいですか?
disp(sol(4).dDT)
27.6483
距離は下限の 25 より大きいですが、下限に近いです。
スケジュールのコストはいくらですか?これを調べるには、exp(sol(4).cost) を計算します。sol(4).cost は総廃棄コストの対数だからです。
disp(exp(sol(4).cost))
2.0883e+07
費用は約2,100万ドルです。
パレート集合内で目標 2 の値が最も低い点を調べます。
この点のインデックスは 3 です。この動作ポイントの金銭的コストははるかに高くなります。
disp(exp(sol(3).cost))
1.4001e+08
金銭的コストは約1億4000万ドルで、これは以前の価値の6倍以上です。しかし、目標 2 が 17 ではなく 10 になることで、廃棄物の埋め立てが 35 (=7*5) 年早まることになります。早めに埋葬する方が安全だと考えられます。
このソリューションの x のスケジュールを表示します。
disp(sol(3).x)
0 0 0 0 0 0 0 0 360 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 240 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 160 0 200 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 240 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 360 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 240 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 360 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 240 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 360 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 240 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 360 0 0 0 0
除去量M(:)と比較して、処分に関するxスケジュールの合計を表示します。
disp([sum(sol(3).x,2),M(:)])
360 360 240 240 360 360 240 240 360 360 240 240 360 360 240 240 360 360 240 240 360 360
繰り返しますが、スケジュールにはすべての削除が考慮されています。
カプセル化施設の稼働時間は何時ですか？
disp(sol(3).ej')
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0
カプセル化は 1 から 18 まで実行されます。
このソリューションの廃棄トンネル間の距離はどれくらいですか?
disp(sol(3).dDT)
38.4942
今回、処分トンネル間の距離は、下限の25メートルと上限の50メートルの中間くらいになります。
まとめ
この例では、問題ベースのアプローチを使用して、非線形、多目的、混合整数最適化問題の定式化を示します。パレート図のデータヒントを使用すると、ソリューションを分析できます。gaplotpareto プロット関数を指定する代わりに、paretoplot 関数を使用してソリューションから同様のプロットを取得できます。
参照
[1] モントネン、アウティ、ティモ・ランタ、マルコ・M・マケラ。対話型多目的最適化による使用済み核燃料の処分スケジュールの計画。アルゴリズム第12巻、第12号、2019年。https://www.mdpi.com/1999-4893/12/12/252 から入手可能。
補助関数
次のコードは、補助関数 max1 を作成します。
function v = max1(A,sij) v = round(max(max(A.*sij - 1))); % "round" ensures integer values end
次のコードは、補助関数 max2 を作成します。
function v = max2(ejOFF) v = round(max((ejOFF').*(1:length(ejOFF)) - 1)); end
参考
gamultiobj | solve | paretoplot | paretosearch
関連するトピック
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)