このページは機械翻訳を使用して翻訳されました。最新版の英語を参照するには、ここをクリックします。
多目的最適化を使用した核燃料処分の計画
この例では、いくつかの整数制約を持つ大規模な非線形多目的問題を定式化して解決する方法を示します。この問題は、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 は、さまざまな適応度関数値を持ついくつかの実行可能解を見つけます。ソリューションに関連付けられた制御変数を見つけるには、データ ヒントを使用します。
データヒントを有効にした後、左上のソリューションをクリックします。データヒントを表示するには、次のコマンドを使用する必要がある場合があります。
ax = gca;
ax.Clipping = "off";
選択されたポイントのインデックスは 8 です。この点に関連付けられた変数は sol(3) にあります。
このソリューションに関連付けられている x 変数を調べます。x(i,j) は、時間 i に削除され、時間 j に処分されることを思い出してください。
disp(sol(8).x)
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 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 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 360 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 0 0 0 0 360 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 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
x スケジュールはすべての削除を考慮しています。
カプセル化施設の稼働時間は何時ですか？
disp(sol(8).ej')
0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0
カプセル化は 14 回目から 18 回目まで実行されます。
廃棄トンネル間の距離はどれくらいですか?
disp(sol(8).dDT)
33.4536
距離は下限の 25 メートルと上限の 50 メートルのほぼ中間です。
スケジュールのコストはいくらですか?これを調べるには、exp(sol(4).cost) を計算します。sol(4).cost は総廃棄コストの対数だからです。
disp(exp(sol(8).cost))
3.1936e+07
費用は約3,200万ドルです。
パレート集合内で目標 2 の値が最も低い点を調べます。
このポイントのインデックスは 10 です。この動作ポイントの金銭的コストははるかに高くなります。
disp(exp(sol(10).cost))
1.2827e+08
金銭的なコストは約1億2800万ドルで、これは従来の4倍以上です。しかし、目標 2 が 17 ではなく 9 になることで、廃棄物の埋め立てが 40 (=8*5) 年早まることになります。早めに埋葬する方が安全だと考えられます。
このソリューションの x のスケジュールを表示します。
disp(sol(10).x)
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 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 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 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 18.4809 37.6345 0 183.8847 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 0 0 0
解はすべて整数値ではありません。除去量 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
繰り返しますが、スケジュールにはすべての削除が考慮されています。
カプセル化施設の稼働時間は何時ですか？
disp(sol(10).ej')
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0
カプセル化は 1 回から 16 回まで実行されます。
このソリューションの廃棄トンネル間の距離はどれくらいですか?
disp(sol(10).dDT)
25
今回は下限距離25メートル。
まとめ
この例では、問題ベースのアプローチを使用して、非線形、多目的、混合整数最適化問題の定式化を示します。パレート図のデータヒントを使用すると、ソリューションを分析できます。gaplotpareto プロット関数を指定する代わりに、paretoplot 関数を使用してソリューションから同様のプロットを取得できます。
参考文献
[1] Montonen, Outi, Timo Ranta, and Marko M. Mäkelä.Planning the Schedule for the Disposal of the Spent Nuclear Fuel with Interactive Multiobjective Optimization.Algorithms Vol. 12, Issue 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
