合金の状態図計算⑤：CALPHAD法による状態図計算【MATLABコード付き】

%% 初期化処理
clc; close all;
clear all;

%% パラメーター設定
dxB = 0.001;
xB = 0.0001 : dxB : 0.9999; % B原子のモル分率
T = 200 : 5 : 1600;         % 計算温度範囲
omega_alpha = 2000;         % α相の原子間相互作用パラメーター
omega_L = 0;                % 液相の原子間相互作用パラメーター
R = 1.987;                  % 気体定数[cal/mol deg]
TA = 900;                   % 原子Aの融点
TB = 1300;                  % 原子Bの融点

% 変数宣言
dG_max = 0;                 % ギブスの自由エネルギー曲線の最大値。可視化範囲決定用。
dG_min = 0;                 % ギブスの自由エネルギー曲線の最小値。可視化範囲決定用。
xmat = [];                  % 最適化結果格納用の変数
xmat_Equilibrium = [];      % 最適化結果格納用の変数（安定平衡線のみ）
xmat_QuasiEquilibrium = []; % 最適化結果格納用の変数（準安定平衡線のみ）

%% 最適化による共通接線探索

% /////////////////////////////////////////
% pflag : 共通接線を探索する自由エネルギー曲線を選択するフラグ
% pflag = 1 : α相と液相の間の共通接線を求める
% pflag = 2 : α相内の共通接線を求める
% pflag = 3 : 液相内の共通接線を求める
% /////////////////////////////////////////

% α相と液相の間の共通接線を求める
pflag = 1;
num_trials = 50; % 異なる初期値で最適化を実行する回数
[xmat, dG_max, dG_min] = ...
    Calphad_optmodule(xB, T, R, TA, TB, omega_alpha, omega_L, xmat, dG_max, dG_min, pflag, num_trials);

% α相内の共通接線を求める
pflag = 2;
num_trials = 50; % 異なる初期値で最適化を実行する回数
[xmat, dG_max, dG_min] = ...
    Calphad_optmodule(xB, T, R, TA, TB, omega_alpha, omega_L, xmat, dG_max, dG_min, pflag, num_trials);

% 液相内の共通接線を求める
pflag = 3;
num_trials = 50; % 異なる初期値で最適化を実行する回数
[xmat, dG_max, dG_min] = ...
    Calphad_optmodule(xB, T, R, TA, TB, omega_alpha, omega_L, xmat, dG_max, dG_min, pflag, num_trials);

%% 可視化

% 可視化用のマトリクスxmat2を構築。
% 解が見つからなかった温度をNanで埋める。
xmat2 = [];
for j = 1 : length(T)
    t = T(j);
    % xmat の中で該当する温度の行を見つける
    rows = xmat(xmat(:, 1) == t, :);
    if isempty(rows)
        % 対応する行がない場合、NaN を使用して温度のみを追加
        xmat2 = [xmat2; t, NaN, NaN, NaN]; % NaNを追加する列数はxmatの列数に合わせる
    else
        % 対応する全ての行を追加
        xmat2 = [xmat2; rows];
    end
end

% 安定平衡線と準安定平衡線の分離計算
for i = 1 : size(xmat2, 1)
    [xmat_Equilibrium, xmat_QuasiEquilibrium] ...
        = SeparateQuasiEquilibrium(xmat2, xB, R, omega_alpha, omega_L, TA, TB, ...
        dG_max, dG_min, i, T, xmat_Equilibrium, xmat_QuasiEquilibrium);
end

% 動画保存
for i = 1 : size(xmat2, 1)
    if i == 1
        figure(2);
        set(gcf, 'Position', [100, 100, 1200, 600]); % ウィンドウサイズを設定
        % plotconfig(xmat, xB, R, omega_alpha, omega_L, TA, TB, dG_max, dG_min, i, T)
        plotconfig(xmat_Equilibrium, xB, R, omega_alpha, omega_L, TA, TB, dG_max, dG_min, i, T)
        filename = ['omega_alpha = ', num2str(omega_alpha, '%.0f'), ', omega_L = ', num2str(omega_L, '%.0f'), '.mp4'];
        v = VideoWriter(filename,'MPEG-4');
        v.FrameRate = 20;
        open(v);
    else
        % plotconfig(xmat, xB, R, omega_alpha, omega_L, TA, TB, dG_max, dG_min, i, T)
        plotconfig(xmat_Equilibrium, xB, R, omega_alpha, omega_L, TA, TB, dG_max, dG_min, i, T)
        frame = getframe(gcf);
        writeVideo(v,frame);
    end
end
close(v);% 動画ファイルを閉じる


%% 以下ローカル関数

% /////////////////////////////////////////
% 可視化設定
% /////////////////////////////////////////
function [] = plotconfig(xmat, xB, R, omega_alpha, omega_L, TA, TB, dG_max, dG_min, i, T)

% 最適化結果から必要な情報を取り出す
Tx = xmat(i, 1);
x1 = xmat(i, 2);
x2 = xmat(i, 3);
pflag = xmat(i, 4);

if pflag == 1 || isnan(pflag)% α相と液相の共通接線の場合、または共通接線を持たない場合。

    % 共通接線の構築
    c0 = (CalcGibssAlpha(omega_alpha, x1, R, Tx) - CalcGibssLiquid(R, TA, x2, TB, omega_L, Tx))/(x1 - x2);
    c1 = CalcGibssAlpha(omega_alpha, x1, R, Tx) - c0 * x1;
    yc = c0 .* xB + c1;

    % 各相のギブスの自由エネルギーの計算
    dG_Alpha = CalcGibssAlpha(omega_alpha, xB, R, Tx);
    dG_L = CalcGibssLiquid(R, TA, xB, TB, omega_L, Tx);

    % 共通接線の描画
    subplot(1,2,1)
    plot(xB, dG_Alpha, xB, dG_L, xB, yc); hold on
    plot(x1, CalcGibssAlpha(omega_alpha, x1, R, Tx), 'ro')
    plot(x2, CalcGibssLiquid(R, TA, x2, TB, omega_L, Tx), 'bo')

elseif pflag == 2 % α相内の共通接線の場合

    % 共通接線の構築
    c0 = (CalcGibssAlpha(omega_alpha, x1, R, Tx) - CalcGibssAlpha(omega_alpha, x2, R, Tx))/(x1 - x2);
    c1 = CalcGibssAlpha(omega_alpha, x1, R, Tx) - c0 * x1;
    yc = c0 .* xB + c1;

    % 各相のギブスの自由エネルギーの計算
    dG_Alpha = CalcGibssAlpha(omega_alpha, xB, R, Tx);
    dG_L = CalcGibssLiquid(R, TA, xB, TB, omega_L, Tx);

    % 共通接線の描画
    subplot(1,2,1)
    plot(xB, dG_Alpha, xB, dG_L, xB, yc); hold on
    plot(x1, CalcGibssAlpha(omega_alpha, x1, R, Tx), 'ro')
    plot(x2, CalcGibssAlpha(omega_alpha, x2, R, Tx), 'bo')

elseif pflag == 3 % 液相内の共通接線の場合

    % 共通接線の構築
    c0 = (CalcGibssLiquid(R, TA, x1, TB, omega_L, Tx)- CalcGibssLiquid(R, TA, x2, TB, omega_L, Tx))/(x1 - x2);
    c1 = CalcGibssLiquid(R, TA, x1, TB, omega_L, Tx) - c0 * x1;
    yc = c0 .* xB + c1;

    % 各相のギブスの自由エネルギーの計算
    dG_Alpha = CalcGibssAlpha(omega_alpha, xB, R, Tx);
    dG_L = CalcGibssLiquid(R, TA, xB, TB, omega_L, Tx);

    % 共通接線の描画
    subplot(1,2,1)
    plot(xB, dG_Alpha, xB, dG_L, xB, yc); hold on
    plot(x1, CalcGibssLiquid(R, TA, x1, TB, omega_L, Tx), 'ro')
    plot(x2, CalcGibssLiquid(R, TA, x2, TB, omega_L, Tx), 'bo')

end

% グラフオプション
title(['temperature = ', num2str(Tx, '%.0f'), ' [K]']);
set(gca, 'FontName', 'Times New Roman', 'FontSize', 12);
%axis equal;
axis tight;
axis on;
fig=gcf;
fig.Color='white';
xlim([0 1])
deltaG = dG_max - dG_min;
ylim([dG_min - deltaG / 10  dG_max + deltaG / 10])
xlabel('x_B')
ylabel('Gibbs free energy')
legend('Galpha','Gliquid','Common Tangent','FontSize', 10)

% 新しいプロットの時、軸設定を保持したまま前のグラフィックスオブジェクトを消去
set(gca,'nextplot','replacechildren');

% 相図の描画

subplot(1,2,2)
red_points = xmat(1 : i, 2);  % 赤い点のx座標
blue_points = xmat(1 : i, 3); % 青い点のx座標
temperatures = xmat(1 : i, 1); % 温度

% 赤い点と青い点をプロット
plot(red_points, temperatures, 'r.'); hold on;
plot(blue_points, temperatures, 'b.');

% グラフオプション
xlim([0 1])
ylim([200 1400])
ylim([min(T) - 50  max(T) + 50])

xlabel('x_B')
ylabel('temperature[K]')
title(['temperature = ', num2str(Tx, '%.0f'), ' [K]']);
set(gca, 'FontName', 'Times New Roman', 'FontSize', 12);

% 新しいプロットの時、軸設定を保持したまま前のグラフィックスオブジェクトを消去
set(gca,'nextplot','replacechildren');

end

% /////////////////////////////////////////
% 目的関数の定義
% /////////////////////////////////////////
function y = obj_fit(x, omega_alpha, R, T, TA, TB,omega_L,pflag)

if pflag == 1

    % α相と液相の間の共通接線を求める目的関数を設定する
    c0 = (CalcGibssAlpha(omega_alpha, x(1), R, T) - CalcGibssLiquid(R, TA, x(2), TB, omega_L, T))/(x(1) - x(2));
    dG_Alpha_prime = CalcGibssAlphaPrime(omega_alpha, x(1), R, T);
    dG_L_prime = CalcGibssLiquidPrime(R, TA, x(2), TB, omega_L, T);
    y = abs(dG_Alpha_prime - c0) + abs(dG_L_prime - c0);  % 目的関数

elseif pflag == 2

    % α相内の共通接線を求める目的関数を設定する
    c0 = (CalcGibssAlpha(omega_alpha, x(1), R, T) - CalcGibssAlpha(omega_alpha, x(2), R, T))/(x(1) - x(2));
    dG_Alpha_prime1 = CalcGibssAlphaPrime(omega_alpha, x(1), R, T);
    dG_Alpha_prime2 = CalcGibssAlphaPrime(omega_alpha, x(2), R, T);
    y = abs(dG_Alpha_prime1 - c0) + abs(dG_Alpha_prime2 - c0);  % 目的関数

    %x(1) と x(2) の差が minDiff 未満の場合にペナルティを適用
    penaltyFactor = CalcPenaltyFactor(x);
    y = y + penaltyFactor;

elseif pflag == 3

    % 液相内の共通接線を求める目的関数を設定する
    c0 = (CalcGibssLiquid(R, TA, x(1), TB, omega_L, T) - CalcGibssLiquid(R, TA, x(2), TB, omega_L, T))/(x(1) - x(2));
    dG_L_prime1 = CalcGibssLiquidPrime(R, TA, x(1), TB, omega_L, T);
    dG_L_prime2 = CalcGibssLiquidPrime(R, TA, x(2), TB, omega_L, T);
    y = abs(dG_L_prime1 - c0) + abs(dG_L_prime2 - c0);  % 目的関数

    %x(1) と x(2) の差が minDiff 未満の場合にペナルティを適用
    penaltyFactor = CalcPenaltyFactor(x);
    y = y + penaltyFactor;

end

end

% /////////////////////////////////////////
% α相のギブスの自由エネルギー
% /////////////////////////////////////////
function dG_alpha = CalcGibssAlpha(omega_alpha, xB, R, T)
dG_alpha = omega_alpha * (1 - xB) .* xB + R * T * ((1 - xB) .* log(1 - xB) + xB .* log(xB));
end

% /////////////////////////////////////////
% 液相のギブスの自由エネルギー
% /////////////////////////////////////////
function dG_L = CalcGibssLiquid(R, Ta, xB, Tb, omega_L, T)
dG_L = R * Ta * (1 - xB) + R * Tb * xB + omega_L * (1 - xB) .* xB - R * T ...
    + R * T * ((1 - xB) .* log(1 - xB) + xB .* log(xB));
end

% /////////////////////////////////////////
% α相のギブスの自由エネルギーの一階微分
% /////////////////////////////////////////
function dG_alpha_p = CalcGibssAlphaPrime(omega_alpha, xB, R, T)
dG_alpha_p = omega_alpha * (1 - 2 .* xB)  + R * T * log(xB./(1 - xB));
end

% /////////////////////////////////////////
% 液相のギブスの自由エネルギーの一階微分
% /////////////////////////////////////////
function dG_L_p = CalcGibssLiquidPrime(R, Ta, xB, Tb, omega_L, T)
dG_L_p = - R * (Ta - Tb) + omega_L * (1 - 2 .* xB)  + R * T * log(xB./(1 - xB));
end

% /////////////////////////////////////////
% 最適化による共通接線の探索
% /////////////////////////////////////////
function [xmat, dG_max, dG_min] =...
    Calphad_optmodule(xB, T, R, TA, TB,omega_alpha, omega_L, xmat, dG_max, dG_min, pflag, num_trials)

% 最適化パラメーター
nvars = 2;           % 変数の数
LB = [0.001 0.001];  % 下限
UB = [0.999 0.999];  % 上限
solutions = zeros(num_trials, nvars);   % 最適解を保存する配列
fvals = zeros(num_trials, 1);           % 目的関数値を保存する配列

% 線形制約（なし）
A = [];
b = [];
Aeq = [];
beq = [];

% 非線形制約（なし）
nonlcon = [];

% 最適化オプションの設定。アルゴリズムは逐次二次計画法（sqp）を選択。
%options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'off');
options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'off', ...
    'TolFun', 1e-6, 'TolX', 1e-6, 'MaxIter', 10000);

% 温度毎に共通接線を探索
for i = 1 : size(T, 2)

    % 進行度の出力
    disp([num2str(T(i)),'[K]/',num2str(T(end)), '[K]を計算中です。']);

    % 可視化範囲決定用に最大値と最小値を求める
    dG_max = max([CalcGibssAlpha(omega_alpha, xB, R, T(i)) ...
        CalcGibssLiquid(R, TA, xB, TB, omega_L, T(i))  dG_max]) ;
    dG_min = min([CalcGibssAlpha(omega_alpha, xB, R, T(i)) ...
        CalcGibssLiquid(R, TA, xB, TB, omega_L, T(i)), dG_min]) ;

    % 目的関関数の設定
    ObjFunc = @(x) obj_fit(x, omega_alpha, R, T(i), TA, TB, omega_L, pflag);

    % 初期値の設定（ラテン超方格）
    [samplingPoints] = CreateInitialPoint(num_trials);

    % 初期値を変更し最適化を実行
    for j = 1 : num_trials

        % 初期値
        % x = [rand(1, 2)]; % 2変数の場合、それぞれの値は 0 から 1 の範囲でランダムに選ばれる
        x = [samplingPoints(j, :)]; % 2変数の場合、それぞれの値は 0 から 1 の範囲でランダムに選ばれる

        % 最適化の実行
        [x, fval] = fmincon(ObjFunc, x, A, b, Aeq, beq, LB, UB, nonlcon, options);

        % 結果の保存
        %fval
        solutions(j, :) = x;
        fvals(j,:) = fval;

    end

    % 目的関数の値が１より大きい解は消去する。
    threshold = 1; % 閾値の設定
    evalmat = [solutions fvals];
    evalmat(evalmat(:, nvars + 1) > threshold, :) = [];
    solutions = evalmat(:, 1 : 2);
    fvals = evalmat(:, 3);

    % 同じ相の自由エネルギー曲線の中で、共通接線の探索を行う場合（pflag=2,3）、
    % 成分が入れ替わっただけの解を見分けるため、昇順に並べ替える。
    % 各行の1列目と2列目を比較し、小さい方を1行目に、大きい方を2行目に移動する。
    if pflag == 2 || pflag == 3
        for l = 1:size(solutions, 1)
            solutions(l, 1:2) = sort(solutions(l, 1:2));
        end
    end

    % 最適化結果の評価
    if size(evalmat, 1) == 0
        disp('失敗。共通接線が見つかりませんでした。');
    else
        disp('成功。共通接線が見つかりました。');

        % クラスタリングで類似した解を一つにまとめる。
        epsilon = 0.01; % 近傍の半径
        MinPts = 1; % クラスタを形成するために必要な最小点数
        idx = dbscan(solutions, epsilon, MinPts);% クラスタリング（DBSCAN）の実行

        % クラスタリング結果の表示
        %disp('クラスタリング結果：');
        %disp(idx);

        % 各クラスタの中心を求める
        unique_clusters = unique(idx);
        cluster_centers = zeros(length(unique_clusters) - 1, size(solutions, 2)); % -1 to exclude noise
        for k = 1 : length(unique_clusters)
            if unique_clusters(k) ~= -1 % ノイズでない場合
                cluster_centers(k, :) = mean(solutions(idx == unique_clusters(k), :), 1);
            end
        end
        disp('クラスタ中心：');
        disp(cluster_centers);

        % 求めた最適値を保存
        Tvec = T(i) * ones(size(cluster_centers,1),1);
        pflagvec = pflag * ones(size(cluster_centers,1),1);
        xtemp = [Tvec cluster_centers pflagvec] ;
        xmat = cat(1, xmat, xtemp);% ここまでの最適値の結果に加える
    end

end

end

% /////////////////////////////////////////
% ペナルティの計算
% /////////////////////////////////////////
function penaltyFactor = CalcPenaltyFactor(x)
minDiff = 0.05; % x(1)とx(2)の最小差異
if abs(x(1) - x(2)) < minDiff
    penaltyFactor = 10^16; % ペナルティの大きさ
else
    penaltyFactor = 0;
end
end

% /////////////////////////////////////////
% 安定平衡線と準安定平衡線の分離計算
% /////////////////////////////////////////
function [xmat_Equilibrium, xmat_QuasiEquilibrium] ...
    = SeparateQuasiEquilibrium(xmat, xB, R, omega_alpha, omega_L, TA, TB, ...
    dG_max, dG_min, i, T, xmat_Equilibrium, xmat_QuasiEquilibrium)

% 最適化結果から必要な情報を取り出す
Tx = xmat(i, 1);
x1 = xmat(i, 2);
x2 = xmat(i, 3);
pflag = xmat(i, 4);

if pflag == 1 || isnan(pflag)% α相と液相の共通接線の場合、または共通接線を持たない場合。

    % 共通接線の構築
    c0 = (CalcGibssAlpha(omega_alpha, x1, R, Tx) - CalcGibssLiquid(R, TA, x2, TB, omega_L, Tx))/(x1 - x2);
    c1 = CalcGibssAlpha(omega_alpha, x1, R, Tx) - c0 * x1;
    yc = c0 .* xB + c1;

    % 各相のギブスの自由エネルギーの計算
    dG_Alpha = CalcGibssAlpha(omega_alpha, xB, R, Tx);
    dG_L = CalcGibssLiquid(R, TA, xB, TB, omega_L, Tx);

elseif pflag == 2 % α相内の共通接線の場合

    % 共通接線の構築
    c0 = (CalcGibssAlpha(omega_alpha, x1, R, Tx) - CalcGibssAlpha(omega_alpha, x2, R, Tx))/(x1 - x2);
    c1 = CalcGibssAlpha(omega_alpha, x1, R, Tx) - c0 * x1;
    yc = c0 .* xB + c1;

    % 各相のギブスの自由エネルギーの計算
    dG_Alpha = CalcGibssAlpha(omega_alpha, xB, R, Tx);
    dG_L = CalcGibssLiquid(R, TA, xB, TB, omega_L, Tx);

elseif pflag == 3 % 液相内の共通接線の場合

    % 共通接線の構築
    c0 = (CalcGibssLiquid(R, TA, x1, TB, omega_L, Tx)- CalcGibssLiquid(R, TA, x2, TB, omega_L, Tx))/(x1 - x2);
    c1 = CalcGibssLiquid(R, TA, x1, TB, omega_L, Tx) - c0 * x1;
    yc = c0 .* xB + c1;

    % 各相のギブスの自由エネルギーの計算
    dG_Alpha = CalcGibssAlpha(omega_alpha, xB, R, Tx);
    dG_L = CalcGibssLiquid(R, TA, xB, TB, omega_L, Tx);

end

% 求めた共通接線が各相のギブスの自由エネルギーよりも小さいかを確認する。
dEmax = max([max(yc - dG_Alpha) max(yc - dG_L)]);

% エネルギーの大小によって配列を振り分ける
if dEmax > 0
    xmat_QuasiEquilibrium(i, :) = xmat(i, :);
    xmat_Equilibrium(i, :) = [Tx, NaN, NaN, NaN];
elseif dEmax <= 0
    xmat_Equilibrium(i, :) = xmat(i, :);
    xmat_QuasiEquilibrium(i, :) = [Tx, NaN, NaN, NaN];
else
    xmat_Equilibrium(i, :) = [Tx, NaN, NaN, NaN];
    xmat_QuasiEquilibrium(i, :) = [Tx, NaN, NaN, NaN];
end

end

% /////////////////////////////////////////
% ラテン超方格法による初期値生成
% /////////////////////////////////////////
function [samplingPoints] = CreateInitialPoint(numPoints)

% 任意の範囲を設定
rangeX = [0.001, 0.999];
rangeY = [0.001, 0.999];

% 2次元のラテン超方格サンプリングを生成
lhsMatrix = lhsdesign(numPoints, 2);

% 範囲に合わせてスケーリング
xPoints = rangeX(1) + (rangeX(2) - rangeX(1)) * lhsMatrix(:, 1);
yPoints = rangeY(1) + (rangeY(2) - rangeY(1)) * lhsMatrix(:, 2);

% XとYが同じになる点を除外し、新しい点を追加
tolerance = 1e-6; % 同一と見なすための許容誤差
i = 1;
while i <= length(xPoints)
    if abs(xPoints(i) - yPoints(i)) < tolerance
        % 新しい点を追加
        newX = rangeX(1) + (rangeX(2) - rangeX(1)) * rand();
        newY = rangeY(1) + (rangeY(2) - rangeY(1)) * rand();
        if abs(newX - newY) >= tolerance
            xPoints(i) = newX;
            yPoints(i) = newY;
        else
            continue; % 新しい点が条件を満たさない場合はスキップ
        end
    end
    i = i + 1;
end

% 点を組み合わせて一つの行列にする
samplingPoints = [xPoints, yPoints];

% % 行列を表示
% disp(samplingPoints);
%
% % 点をプロット
% scatter(samplingPoints(:, 1), samplingPoints(:, 2));
% title(['', num2str(numPoints), '点のラテン超方格サンプリング（X ≠ Y）']);
% xlabel('X軸');
% ylabel('Y軸');

end