#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ 一ダイオードモデル（SDM）を用いた太陽電池IV特性のフィッティングツール。 このスクリプトは、測定された電流-電圧（IV）特性データに対し、 一ダイオードモデルに基づいてフィッティングを行い、モデルパラメータを抽出します。 初期パラメータの推定、最適化、結果の保存、およびIV曲線のプロット機能を提供します。 関連リンク: :doc:`pvfit_usage` """ import argparse import csv import math import os import sys import traceback from pathlib import Path import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize, brentq # --- 定数 --- KB = 1.380649e-23 E_CHARGE = 1.602176634e-19 PARAM_NAMES = ["I0", "ndiode", "IPV", "Rs", "Rsh"] EPS_I = 1.0e-15 EPS_R = 1.0e-6 # --- データ入出力 --- def read_data(infile, xmin=None, xmax=None): """ 入力CSVファイルからIV特性データを読み込み、電圧(V)の昇順でソートして返します。 データはCSVファイルから読み込まれ、"DataValue"セクションに続く数値がIVデータとして解釈されます。 オプションで電圧範囲を指定してデータをフィルタリングできます。 データが存在しない場合や無効なデータが多い場合はエラーを発生させます。 :param infile: 読み込むCSVファイルのパス。 :type infile: str :param xmin: プロットおよびフィッティングに使用する電圧範囲の最小値 (V)。Noneの場合、制限なし。 :type xmin: float or None :param xmax: プロットおよびフィッティングに使用する電圧範囲の最大値 (V)。Noneの場合、制限なし。 :type xmax: float or None :returns: - V (np.ndarray): 読み込まれた電圧データ（ソート済み）。 - I (np.ndarray): 読み込まれた電流データ（ソート済み）。 - inf (dict): ファイル名と記録時間を含む情報辞書。 :rtype: tuple[np.ndarray, np.ndarray, dict] :raises ValueError: 有効なデータがファイル内に見つからない場合。 """ data_started = False raw_x, raw_y = [], [] inf = {"RecordTime": "Unknown", "FileName": infile} with open(infile, "r", encoding="utf-8-sig", newline="") as f: reader = csv.reader(f) for row in reader: if not row: continue cols = [c.strip() for c in row] if len(cols) < 2: continue if cols[0] == "MetaData" and cols[1] == "TestRecord.RecordTime": if len(cols) >= 3: inf["RecordTime"] = cols[2] if cols[0] == "DataName": data_started = True continue if data_started and cols[0] == "DataValue": try: raw_x.append(float(cols[1])) raw_y.append(float(cols[2])) except (ValueError, IndexError): continue if not raw_x: raise ValueError(f"No valid data found in {infile}") # 配列化 x_arr, y_arr = np.array(raw_x), np.array(raw_y) # --- 重要: 電圧 V でソートを行う --- # Double Sweep（往復測定）データでも fill_between が正しく動くようにするため sort_idx = np.argsort(x_arr) x_arr = x_arr[sort_idx] y_arr = y_arr[sort_idx] # --------------------------------- print(f"File Name: {inf['FileName']}") print(f"Record Time: {inf['RecordTime']}") print(f"Voltage range: {x_arr.min():.4f} to {x_arr.max():.4f} V") mask = np.ones(len(x_arr), dtype=bool) if xmin is not None: mask &= (x_arr >= xmin) if xmax is not None: mask &= (x_arr <= xmax) return x_arr[mask], y_arr[mask], inf def load_param_csv(csv_path): """ 指定されたCSVファイルからモデルパラメータと固定するパラメータのセットを読み込みます。 CSVファイルには、"varname", "value", "optid" (0:固定, 1:自由) の列が必要です。 `PARAM_NAMES` に含まれるパラメータのみがフィッティングの対象として扱われ、 `optid` が0のパラメータは固定としてマークされます。 それ以外のパラメータは、フィッティング対象外の補助情報として読み込まれます。 :param csv_path: パラメータを読み込むCSVファイルのパス。 :type csv_path: str :returns: - params (dict): 読み込まれたパラメータ名とその値の辞書。 - fix_set (set): 固定としてマークされたパラメータ名のセット。 :rtype: tuple[dict, set] """ params, fix_set = {}, set() if os.path.exists(csv_path): with open(csv_path, "r", encoding="utf-8") as f: reader = csv.DictReader(f) for row in reader: name = row["varname"] try: val = float(row["value"]) if name in PARAM_NAMES: params[name] = val if int(row.get("optid", 1)) == 0: fix_set.add(name) else: params[name] = val # フィッティングパラメータ以外の値も読み込む except (ValueError, TypeError): continue return params, fix_set def save_param_csv(csv_path, params, fix_set, errors=None, rss=None): """ モデルパラメータ、固定設定、推定誤差、およびRSSをCSVファイルに保存します。 結果として得られたフィッティングパラメータ、その固定/自由設定、 そして推定された誤差（利用可能な場合）を整理してCSV形式で出力します。 推定された開回路電圧 (VOC_est) と短絡電流 (ISC_est) も含められます。 :param csv_path: パラメータを保存するCSVファイルのパス。 :type csv_path: str :param params: 保存するパラメータ名とその値の辞書。 :type params: dict :param fix_set: 固定されたパラメータ名のセット。 :type fix_set: set :param errors: 各パラメータの推定誤差の辞書。デフォルトはNone。 :type errors: dict or None :param rss: 残差平方和 (Residual Sum of Squares) の値。デフォルトはNone。 :type rss: float or None :returns: なし :rtype: None """ if errors is None: errors = {} with open(csv_path, "w", encoding="utf-8", newline="") as f: writer = csv.writer(f) writer.writerow(["varname", "value", "optid", "error"]) for name in PARAM_NAMES: err = errors.get(name, 0) if errors else 0 writer.writerow([name, params[name], 0 if name in fix_set else 1, err]) writer.writerow(["VOC_est", params.get("VOC", 0), 0, 0]) writer.writerow(["ISC_est", params.get("ISC", 0), 0, 0]) if rss is not None: writer.writerow(["RSS_logy", rss, 0, 0]) # --- 物理モデル・計算 --- def solve_root_poly(x, y, target_y=0, order=3): """ 多項式近似を用いて、与えられたデータセット(x, y)からyが特定の`target_y`になるxを推定します。 与えられたIVデータに対して、n次多項式で近似し、`y - target_y = 0`となる根を探します。 複数の実数根が見つかった場合、データの平均x値に最も近い根を選択します。 実数根が見つからない場合は、線形補間を使用します。 :param x: x軸のデータポイント（電圧）。 :type x: np.ndarray :param y: y軸のデータポイント（電流）。 :type y: np.ndarray :param target_y: 求めるyの値。デフォルトは0（開回路電圧VOCの推定に相当）。 :type target_y: float :param order: 多項式近似の次数。デフォルトは3。 :type order: int :returns: `target_y` に対応する推定されたx値。 :rtype: float """ if len(x) < order + 1: order = max(1, len(x) - 1) z = np.polyfit(x, y - target_y, order) roots = np.roots(z) real_roots = roots[np.isreal(roots)].real if len(real_roots) == 0: return np.interp(target_y, y, x) return real_roots[np.argmin(np.abs(real_roots - np.mean(x)))] def model(V, I0, ndiode, IPV, Rs, Rsh, temperature=300.0): """ 一ダイオードモデル（Single-Diode Model, SDM）に基づいてIV特性の電流を計算します。 このモデルは以下の数式で表されます： I = IPV - I0 * (exp((V + I*Rs) / (ndiode * k * T / q)) - 1) - (V + I*Rs) / Rsh ここで、Iは全電流、IPVは光電流、I0は飽和電流、ndiodeはダイオード理想因子、 Rsは直列抵抗、Rshはシャント抵抗、kはボルツマン定数、Tは温度、qは素電荷です。 この式はIについて解析的に解けないため、`scipy.optimize.brentq` を用いて数値的に解を求めます。 :param V: 電圧のNumpy配列または単一値。 :type V: np.ndarray or float :param I0: 飽和電流 (A)。 :type I0: float :param ndiode: ダイオード理想因子。 :type ndiode: float :param IPV: 光電流 (A)。 :type IPV: float :param Rs: 直列抵抗 (Ω)。 :type Rs: float :param Rsh: シャント抵抗 (Ω)。 :type Rsh: float :param temperature: セル温度 (K)。デフォルトは300.0。 :type temperature: float :returns: 入力電圧Vに対応する計算された電流 (A) のNumpy配列。 :rtype: np.ndarray """ V = np.atleast_1d(V) vt = ndiode * KB * temperature / E_CHARGE def resid(Idark, v_val): vd = v_val - Rs * Idark arg = np.clip(vd / vt, -700, 700) # オーバーフロー防止 return I0 * (np.exp(arg) - 1.0) + vd / Rsh - Idark results = [] for v_val in V: try: # 探索範囲を適切に設定 limit = abs(v_val / max(Rsh, EPS_R)) + abs(I0) * 1e5 + abs(IPV) + 1.0 id_val = brentq(resid, -limit, limit, args=(v_val,)) results.append(id_val - IPV) except (ValueError, RuntimeError): # brentqが根を見つけられない場合 # フォールバックとして、Rsが非常に大きい場合の近似を返す # または、単純な抵抗モデルに近い値を返す results.append(v_val / max(Rsh, EPS_R) - IPV) return np.array(results) def objective(p_free, V, I_meas, temp, base_params, fix_set): """ 最適化アルゴリズムのための目的関数を計算します。 この関数は、自由パラメータ`p_free`を更新し、それらの値を用いてモデルを呼び出し、 測定された電流と計算された電流の対数絶対値の差の二乗和（残差平方和）を返します。 `ndiode`を除くパラメータは対数スケールで最適化されます。 `EPS_I` は、対数計算でのゼロ除算や数値不安定性を防ぐために使用されます。 :param p_free: 最適化中の自由パラメータの配列。 :type p_free: np.ndarray :param V: 測定された電圧データ。 :type V: np.ndarray :param I_meas: 測定された電流データ。 :type I_meas: np.ndarray :param temp: セル温度 (K)。 :type temp: float :param base_params: 固定パラメータを含む、全てのパラメータの初期値または現在の値の辞書。 :type base_params: dict :param fix_set: 固定されたパラメータ名のセット。 :type fix_set: set :returns: 測定値とモデル予測値の対数絶対値の差の二乗和。 :rtype: float """ curr = dict(base_params); idx = 0 for name in PARAM_NAMES: if name not in fix_set: val = p_free[idx] curr[name] = 10**val if name != "ndiode" else val # ndiode以外はlog10スケールで最適化 idx += 1 I_calc = model(V, **{k: curr[k] for k in PARAM_NAMES}, temperature=temp) # 微小値EPS_Iを加えてログ計算でのゼロ除算を防ぐ return np.sum((np.log10(np.abs(I_meas) + EPS_I) - np.log10(np.abs(I_calc) + EPS_I))**2) # --- 誤差推定関連 --- def get_jacobian(p_free, V, temp, base_params, fix_set): """ モデルの対数電流に対するヤコビ行列を数値的に計算します。 ヤコビ行列は、各自由パラメータに関するモデル出力（対数電流）の偏微分を近似したものです。 有限差分法を用いて計算され、誤差推定に使用されます。 `EPS_I` は、対数計算でのゼロ除算や数値不安定性を防ぐために使用されます。 :param p_free: 現在の自由パラメータの配列。 :type p_free: np.ndarray :param V: 測定された電圧データ。 :type V: np.ndarray :param temp: セル温度 (K)。 :type temp: float :param base_params: 固定パラメータを含む、全てのパラメータの初期値または現在の値の辞書。 :type base_params: dict :param fix_set: 固定されたパラメータ名のセット。 :type fix_set: set :returns: モデルの対数電流に対するヤコビ行列。 :rtype: np.ndarray """ eps = 1e-5 # 微小な摂動値 n_params, n_points = len(p_free), len(V) J = np.zeros((n_points, n_params)) def get_log_model(p): """与えられたパラメータでモデル電流を計算し、その対数絶対値を返すヘルパー関数""" curr = dict(base_params); idx = 0 for name in PARAM_NAMES: if name not in fix_set: val = p[idx] curr[name] = 10**val if name != "ndiode" else val idx += 1 I_c = model(V, **{k: curr[k] for k in PARAM_NAMES}, temperature=temp) return np.log10(np.abs(I_c) + EPS_I) f0 = get_log_model(p_free) for i in range(n_params): p_eps = np.copy(p_free); p_eps[i] += eps # 各パラメータを微小に摂動 J[:, i] = (get_log_model(p_eps) - f0) / eps # 有限差分で偏微分を計算 return J def estimate_errors(res, V, I_meas, temp, base_params, fix_set): """ 非線形最小二乗法で得られた結果に基づき、パラメータの誤差とモデル曲線の標準偏差を推定します。 フィッティング結果（`res`）とヤコビ行列を用いて、線形近似に基づいた共分散行列を計算します。 この共分散行列から、各パラメータの標準誤差と、モデル曲線自体の信頼性区間を推定します。 `np.linalg.pinv` を使用することで、数値的に不安定な場合でも計算を継続します。 パラメータの誤差は、`ndiode`は線形スケール、それ以外は元のスケール（対数スケールから変換）で報告されます。 :param res: `scipy.optimize.minimize` の結果オブジェクト。最適化されたパラメータ (`res.x`) と目的関数値 (`res.fun`) を含む。 :type res: scipy.optimize.OptimizeResult :param V: 測定された電圧データ。 :type V: np.ndarray :param I_meas: 測定された電流データ。 :type I_meas: np.ndarray :param temp: セル温度 (K)。 :type temp: float :param base_params: 固定パラメータを含む、全てのパラメータの初期値または現在の値の辞書。 :type base_params: dict :param fix_set: 固定されたパラメータ名のセット。 :type fix_set: set :returns: - errors_dict (dict): 各自由パラメータの推定誤差の辞書。 - sigma_log_I (np.ndarray): モデル曲線の対数スケールでの標準偏差。 :rtype: tuple[dict, np.ndarray] """ p_opt = res.x n, p = len(I_meas), len(p_opt) if n <= p: # データポイント数がパラメータ数以下の場合、誤差推定は不可能 return {}, np.zeros(len(V)) J = get_jacobian(p_opt, V, temp, base_params, fix_set) sigma2_hat = res.fun / (n - p) # 残差の分散を推定 try: # inv ではなく pinv (擬似逆行列) を使うことで、数値的に不安定な場合でも計算を継続 # J.T @ J は Hessian の近似（Gauss-Newton法における） cov = sigma2_hat * np.linalg.pinv(J.T @ J) p_errors_log = np.sqrt(np.maximum(np.diag(cov), 0)) # パラメータの対数スケールでの標準誤差 sigma_log_I = np.sqrt(np.maximum(np.diag(J @ cov @ J.T), 0)) # モデル曲線の対数スケールでの標準偏差 errors_dict, idx = {}, 0 for name in PARAM_NAMES: if name not in fix_set: if name == "ndiode": errors_dict[name] = p_errors_log[idx] # ndiodeは線形スケール else: # log10スケールの誤差から元のスケールでの誤差に変換 # d(10^x)/dx = 10^x * ln(10) を利用 errors_dict[name] = (10**p_opt[idx]) * np.log(10) * p_errors_log[idx] idx += 1 return errors_dict, sigma_log_I except Exception as e: print(f"Warning: Error estimation failed ({e})") return {}, np.zeros(len(V)) # --- 描画ロジック --- def plot_iv(V, I_meas, I_final, axes, sigma_log=None, label="Model", title=""): """ 測定IV特性とモデルIV特性を線形スケールおよび対数スケールでプロットします。 線形スケールのプロットと対数スケールのプロットの2つのサブプロット（axes）に、 測定データ、フィッティングされたモデル曲線、および（指定された場合）信頼性区間を描画します。 対数スケールプロットでは電流の絶対値が使用され、`EPS_I` を加えることでゼロ値のプロ処理を避けます。 :param V: 電圧データ。 :type V: np.ndarray :param I_meas: 測定された電流データ。 :type I_meas: np.ndarray :param I_final: モデルによって計算された最終電流データ。 :type I_final: np.ndarray :param axes: プロットに使用するmatplotlibの2つのサブプロット軸 (ax_lin, ax_log)。 :type axes: tuple[matplotlib.axes.Axes, matplotlib.axes.Axes] :param sigma_log: モデル曲線の対数スケールでの標準偏差。提供された場合、信頼性区間が描画されます。 :type sigma_log: np.ndarray or None :param label: モデル曲線の凡例ラベル。デフォルトは"Model"。 :type label: str :param title: プロットのタイトル。デフォルトは空文字列。 :type title: str :returns: なし :rtype: None """ ax_lin, ax_log = axes ax_lin.plot(V, I_meas, "o", color="blue", markersize=4, alpha=0.3, label="Data") ax_lin.plot(V, I_final, "-", color="orange", lw=2, label=label) # 信頼性区間の塗りつぶし (水色領域) if sigma_log is not None: upper = 10**(np.log10(np.abs(I_final) + EPS_I) + sigma_log) lower = 10**(np.log10(np.abs(I_final) + EPS_I) - sigma_log) # np.sign(I_final) を乗算することで、電流の符号を維持しつつ信頼性区間をプロット ax_lin.fill_between(V, np.sign(I_final)*lower, np.sign(I_final)*upper, color='lightblue', alpha=0.5) # 対数プロットは絶対値でプロットされるため、符号は不要 ax_log.fill_between(V, lower, upper, color='lightblue', alpha=0.5) ax_lin.axhline(0, color='k', lw=0.5); ax_lin.axvline(0, color='k', lw=0.5) ax_lin.set_title(title); ax_lin.grid(True); ax_lin.legend() # 対数プロットは絶対値で描画されるため、EPS_Iを加えてゼロを防ぐ ax_log.semilogy(V, np.abs(I_meas) + EPS_I, "o", color="blue", markersize=4, alpha=0.3) ax_log.semilogy(V, np.abs(I_final) + EPS_I, "-", color="orange", lw=2) ax_log.grid(True, which="both") # --- モード別実行 --- def run_fit_mode(V, I_meas, init_params, args, csv_path, fix_set): """ フィッティングモードを実行し、IV特性データに一ダイオードモデルを適合させます。 初期パラメータから最適化を開始し、`scipy.optimize.minimize` を用いて目的関数を最小化します。 フィッティングの進行状況は、指定された間隔でプロットとコンソール出力で表示されます。 最終的なフィッティングパラメータ、その誤差、および推定されたVOCとISCをCSVファイルに保存します。 `ndiode`を除くパラメータは対数スケールに変換されて最適化されます。 :param V: 測定された電圧データ。 :type V: np.ndarray :param I_meas: 測定された電流データ。 :type I_meas: np.ndarray :param init_params: フィッティングの初期パラメータの辞書。 :type init_params: dict :param args: コマンドライン引数を格納したargparse.Namespaceオブジェクト。 :type args: argparse.Namespace :param csv_path: 結果のパラメータを保存するCSVファイルのパス。 :type csv_path: str :param fix_set: 固定するパラメータ名のセット。 :type fix_set: set :returns: - fit_params (dict): フィッティング後の最適化されたパラメータの辞書。 - sigma_log_curve (np.ndarray): 最適化されたモデル曲線の対数スケールでの標準偏差。 :rtype: tuple[dict, np.ndarray] """ # ndiode以外はlog10スケールで最適化するため変換 p0_free = [init_params[p] if p == "ndiode" else np.log10(max(init_params[p], EPS_R if "R" in p else EPS_I)) for p in PARAM_NAMES if p not in fix_set] plt.ion(); fig, axes = plt.subplots(1, 2, figsize=(11, 5)); iters = [0] def callback(xk): """最適化中の進捗をプロット・表示するためのコールバック関数""" iters[0] += 1 if iters[0] % args.ninterval_plot == 0: curr = dict(init_params); idx = 0 for name in PARAM_NAMES: if name not in fix_set: curr[name] = 10**xk[idx] if name != "ndiode" else xk[idx]; idx += 1 I_c = model(V, **{k: curr[k] for k in PARAM_NAMES}, temperature=args.temperature) rss = objective(xk, V, I_meas, args.temperature, init_params, fix_set) for ax in axes: ax.cla() # 既存のプロットをクリア plot_iv(V, I_meas, I_c, axes, label=f"Iter {iters[0]}", title=f"RSS: {rss:.2e}") plt.pause(0.01) # 短時間一時停止してGUIを更新 if iters[0] % args.ninterval_print == 0: p_str = " | ".join([f"{n}: {curr[n]:.2e}" + ("(F)" if n in fix_set else "") for n in PARAM_NAMES]) print(f"Iter {iters[0]:4d} | RSS: {rss:.4e} | {p_str}") res = minimize(objective, p0_free, args=(V, I_meas, args.temperature, init_params, fix_set), method=args.method, callback=callback) plt.ioff(); plt.close(fig) # インタラクティブモードを終了し、フィッティング中のプロットウィンドウを閉じる fit_params = dict(init_params); idx = 0 for name in PARAM_NAMES: if name not in fix_set: val = res.x[idx]; fit_params[name] = 10**val if name != "ndiode" else val; idx += 1 errors, sigma_log_curve = estimate_errors(res, V, I_meas, args.temperature, init_params, fix_set) # 推定VOCとISCを計算 # ISCはV=0での電流値 fit_params["ISC"] = model(0, **{k: fit_params[k] for k in PARAM_NAMES}, temperature=args.temperature)[0] # VOCはI=0となる電圧値（多項式根探索で推定） v_f = np.linspace(min(V.min(), 0), max(V.max(), 1.2), 2000) # VOC推定のための広範囲電圧データ fit_params["VOC"] = solve_root_poly(v_f, model(v_f, **{k: fit_params[k] for k in PARAM_NAMES}, temperature=args.temperature)) save_param_csv(csv_path, fit_params, fix_set, errors=errors, rss=res.fun) print(f"\nFinal Results (RSS: {res.fun:.4e})") for p in PARAM_NAMES: err_str = f" +/- {errors[p]:.2e}" if p in errors else " (Fixed)" print(f" {p:8s}: {fit_params[p]:.6e}{err_str}") return fit_params, sigma_log_curve def get_initial_params(V, I_meas, args, csv_path): """ フィッティングのための初期パラメータを生成または既存のCSVファイルから読み込みます。 `--mode init` が指定されている場合、IV曲線の特徴（V=0近傍の電流からIPV、Vが負の領域の傾きからRsh） から初期パラメータを推定します。 それ以外の場合、既存のパラメータCSVファイルから読み込みます。 CSVファイルが存在しない、または空の場合は、自動的に初期パラメータが推定されます。 コマンドライン引数で明示的に指定されたパラメータは、常に読み込み値や推定値を上書きします。 :param V: 測定された電圧データ。 :type V: np.ndarray :param I_meas: 測定された電流データ。 :type I_meas: np.ndarray :param args: コマンドライン引数を格納したargparse.Namespaceオブジェクト。 :type args: argparse.Namespace :param csv_path: 初期パラメータを読み書きするCSVファイルのパス。 :type csv_path: str :returns: - params (dict): 初期パラメータ名とその値の辞書。 - csv_fix (set): CSVファイルで固定として指定されたパラメータ名のセット。 :rtype: tuple[dict, set] """ if args.mode == "init": # 初期パラメータの推定ロジック # V=0近傍の電流からIPVを推定 uV, idxs = np.unique(V, return_index=True); uI = I_meas[idxs] # 重複電圧を除去 idx_v0 = np.argmin(np.abs(uV)) # V=0に最も近いインデックス # V=0近傍のIV曲線から、直線近似/多項式近似でISCを推定（IPVに相当） ipv_est = max(-np.polyval(np.polyfit(uV[max(0,idx_v0-3):idx_v0+4], uI[max(0,idx_v0-3):idx_v0+4], 3), 0), EPS_I) # Vが負の領域の傾きからRshを推定 rsh_est = 1.0 / max(np.gradient(uI, uV)[uV < 0][0], 1e-12) if np.any(uV < 0) else 1e6 params = {"I0": 1e-10, "ndiode": 1.5, "IPV": ipv_est, "Rs": 10.0, "Rsh": rsh_est, "VOC": 0.5, "ISC": -ipv_est} csv_fix = set() else: params, csv_fix = load_param_csv(csv_path) if not params: # CSVファイルが存在しない、または空の場合、initモードで初期化 print(f"Warning: Parameter CSV '{csv_path}' not found or empty. Generating initial parameters.") temp_args = argparse.Namespace(**{**vars(args), 'mode': 'init'}) # 一時的にinitモードに設定 params, _ = get_initial_params(V, I_meas, temp_args, csv_path) # コマンドライン引数で指定されたパラメータで上書き for p in PARAM_NAMES: val = getattr(args, p); if val is not None: params[p] = val return params, csv_fix # --- メイン --- def main(): """ スクリプトの主要な実行フローを制御します。 コマンドライン引数を解析し、指定されたモード（init, fit, sim）に基づいて、 データの読み込み、初期パラメータの取得、フィッティング実行、シミュレーション、 そして結果の保存とプロットを行います。 エラー発生時にはトレースバックを出力し、プログラムを終了します。 """ parser = argparse.ArgumentParser(description="一ダイオードモデルによる太陽電池IV特性フィッティング") parser.add_argument("--mode", choices=["init", "fit", "sim"], default="fit", help="実行モード: init (初期パラメータ生成), fit (フィッティング), sim (シミュレーション)") parser.add_argument("--method", default="Nelder-Mead", help="scipy.optimize.minimizeで使う最適化アルゴリズム (例: Nelder-Mead, L-BFGS-B)") parser.add_argument("--infile", default="input.csv", help="入力IV特性CSVファイルのパス") parser.add_argument("--temperature", type=float, default=300.0, help="計算に使用するセル温度 (K)") parser.add_argument("--xmin", type=float, help="プロットおよびフィッティングに使用する電圧範囲の最小値 (V)") parser.add_argument("--xmax", type=float, help="プロットおよびフィッティングに使用する電圧範囲の最大値 (V)") for p in PARAM_NAMES: # 各モデルパラメータに対する引数を動的に追加 parser.add_argument(f"--{p}", type=float, help=f"パラメータ {p} の値を指定 (既存設定を上書き)") parser.add_argument("--fix", nargs="*", default=[], help="フィッティング中に固定するパラメータの名前 (例: I0 ndiode)") parser.add_argument("--nlsq_points", type=int, default=7, # 現在未使用だが、将来の拡張用 help="非線形最小二乗近似のためのデータポイント数") parser.add_argument("--ninterval_print", type=int, default=10, help="フィッティング中のコンソール出力間隔 (イテレーション数)") parser.add_argument("--ninterval_plot", type=int, default=10, help="フィッティング中のプロット更新間隔 (イテレーション数)") args = parser.parse_args() V, I_meas, inf = read_data(args.infile, xmin=args.xmin, xmax=args.xmax) csv_path = Path(args.infile).stem + "-parameters.csv" params, csv_fix = get_initial_params(V, I_meas, args, csv_path) fix_set = set(args.fix) | csv_fix # コマンドライン引数とCSV設定の両方で固定されたパラメータを結合 sigma_log = None if args.mode == "init": save_param_csv(csv_path, params, fix_set); print("\n--- Initial Parameters Saved ---") elif args.mode == "fit": params, sigma_log = run_fit_mode(V, I_meas, params, args, csv_path, fix_set) elif args.mode == "sim": print("\n--- Simulation Mode ---") for p in PARAM_NAMES: print(f" {p:8s}: {params[p]:.6e}") # 最終結果のプロット I_final = model(V, **{k: params[k] for k in PARAM_NAMES}, temperature=args.temperature) fig, axes = plt.subplots(1, 2, figsize=(11, 5)) plot_iv(V, I_meas, I_final, axes, sigma_log=sigma_log, label=f"Final ({args.mode})", title=inf["FileName"]) plt.tight_layout(); plt.show() if __name__ == "__main__": try: main() except Exception: traceback.print_exc(); sys.exit(1)