pvfit_simple.py ダウンロード/コピー

pvfit_simple.py をダウンロード

pvfit_simple.py
pvfit_simple.py
  1#!/usr/bin/env python3
  2# -*- coding: utf-8 -*-
  3"""
  4一ダイオードモデル（SDM）を用いた太陽電池IV特性のフィッティングツール。
  5
  6このスクリプトは、測定された電流-電圧（IV）特性データに対し、
  7一ダイオードモデルに基づいてフィッティングを行い、モデルパラメータを抽出します。
  8初期パラメータの推定、最適化、結果の保存、およびIV曲線のプロット機能を提供します。
  9
 10関連リンク:
 11:doc:`pvfit_usage`
 12"""
 13
 14import argparse
 15import csv
 16import math
 17import os
 18import sys
 19import traceback
 20from pathlib import Path
 21
 22import numpy as np
 23import matplotlib.pyplot as plt
 24from scipy.optimize import minimize, brentq
 25
 26# --- 定数 ---
 27KB = 1.380649e-23
 28E_CHARGE = 1.602176634e-19
 29PARAM_NAMES = ["I0", "ndiode", "IPV", "Rs", "Rsh"]
 30EPS_I = 1.0e-15
 31EPS_R = 1.0e-6
 32
 33# --- データ入出力 ---
 34
 35def read_data(infile, xmin=None, xmax=None):
 36    """
 37    入力CSVファイルからIV特性データを読み込み、電圧(V)の昇順でソートして返します。
 38
 39    データはCSVファイルから読み込まれ、"DataValue"セクションに続く数値がIVデータとして解釈されます。
 40    オプションで電圧範囲を指定してデータをフィルタリングできます。
 41    データが存在しない場合や無効なデータが多い場合はエラーを発生させます。
 42
 43    :param infile: 読み込むCSVファイルのパス。
 44    :type infile: str
 45    :param xmin: プロットおよびフィッティングに使用する電圧範囲の最小値 (V)。Noneの場合、制限なし。
 46    :type xmin: float or None
 47    :param xmax: プロットおよびフィッティングに使用する電圧範囲の最大値 (V)。Noneの場合、制限なし。
 48    :type xmax: float or None
 49    :returns:
 50        - V (np.ndarray): 読み込まれた電圧データ（ソート済み）。
 51        - I (np.ndarray): 読み込まれた電流データ（ソート済み）。
 52        - inf (dict): ファイル名と記録時間を含む情報辞書。
 53    :rtype: tuple[np.ndarray, np.ndarray, dict]
 54    :raises ValueError: 有効なデータがファイル内に見つからない場合。
 55    """
 56    data_started = False
 57    raw_x, raw_y = [], []
 58    inf = {"RecordTime": "Unknown", "FileName": infile}
 59    
 60    with open(infile, "r", encoding="utf-8-sig", newline="") as f:
 61        reader = csv.reader(f)
 62        for row in reader:
 63            if not row: continue
 64            cols = [c.strip() for c in row]
 65            if len(cols) < 2: continue
 66            if cols[0] == "MetaData" and cols[1] == "TestRecord.RecordTime":
 67                if len(cols) >= 3: inf["RecordTime"] = cols[2]
 68            if cols[0] == "DataName":
 69                data_started = True
 70                continue
 71            if data_started and cols[0] == "DataValue":
 72                try:
 73                    raw_x.append(float(cols[1]))
 74                    raw_y.append(float(cols[2]))
 75                except (ValueError, IndexError): continue
 76    
 77    if not raw_x: raise ValueError(f"No valid data found in {infile}")
 78    
 79    # 配列化
 80    x_arr, y_arr = np.array(raw_x), np.array(raw_y)
 81    
 82    # --- 重要: 電圧 V でソートを行う ---
 83    # Double Sweep（往復測定）データでも fill_between が正しく動くようにするため
 84    sort_idx = np.argsort(x_arr)
 85    x_arr = x_arr[sort_idx]
 86    y_arr = y_arr[sort_idx]
 87    # ---------------------------------
 88
 89    print(f"File Name: {inf['FileName']}")
 90    print(f"Record Time: {inf['RecordTime']}")
 91    print(f"Voltage range: {x_arr.min():.4f} to {x_arr.max():.4f} V")
 92
 93    mask = np.ones(len(x_arr), dtype=bool)
 94    if xmin is not None: mask &= (x_arr >= xmin)
 95    if xmax is not None: mask &= (x_arr <= xmax)
 96    
 97    return x_arr[mask], y_arr[mask], inf
 98
 99
100def load_param_csv(csv_path):
101    """
102    指定されたCSVファイルからモデルパラメータと固定するパラメータのセットを読み込みます。
103
104    CSVファイルには、"varname", "value", "optid" (0:固定, 1:自由) の列が必要です。
105    `PARAM_NAMES` に含まれるパラメータのみがフィッティングの対象として扱われ、
106    `optid` が0のパラメータは固定としてマークされます。
107    それ以外のパラメータは、フィッティング対象外の補助情報として読み込まれます。
108
109    :param csv_path: パラメータを読み込むCSVファイルのパス。
110    :type csv_path: str
111    :returns:
112        - params (dict): 読み込まれたパラメータ名とその値の辞書。
113        - fix_set (set): 固定としてマークされたパラメータ名のセット。
114    :rtype: tuple[dict, set]
115    """
116    params, fix_set = {}, set()
117    if os.path.exists(csv_path):
118        with open(csv_path, "r", encoding="utf-8") as f:
119            reader = csv.DictReader(f)
120            for row in reader:
121                name = row["varname"]
122                try:
123                    val = float(row["value"])
124                    if name in PARAM_NAMES:
125                        params[name] = val
126                        if int(row.get("optid", 1)) == 0: fix_set.add(name)
127                    else: params[name] = val # フィッティングパラメータ以外の値も読み込む
128                except (ValueError, TypeError): continue
129    return params, fix_set
130
131def save_param_csv(csv_path, params, fix_set, errors=None, rss=None):
132    """
133    モデルパラメータ、固定設定、推定誤差、およびRSSをCSVファイルに保存します。
134
135    結果として得られたフィッティングパラメータ、その固定/自由設定、
136    そして推定された誤差（利用可能な場合）を整理してCSV形式で出力します。
137    推定された開回路電圧 (VOC_est) と短絡電流 (ISC_est) も含められます。
138
139    :param csv_path: パラメータを保存するCSVファイルのパス。
140    :type csv_path: str
141    :param params: 保存するパラメータ名とその値の辞書。
142    :type params: dict
143    :param fix_set: 固定されたパラメータ名のセット。
144    :type fix_set: set
145    :param errors: 各パラメータの推定誤差の辞書。デフォルトはNone。
146    :type errors: dict or None
147    :param rss: 残差平方和 (Residual Sum of Squares) の値。デフォルトはNone。
148    :type rss: float or None
149    :returns: なし
150    :rtype: None
151    """
152    if errors is None: errors = {}
153    with open(csv_path, "w", encoding="utf-8", newline="") as f:
154        writer = csv.writer(f)
155        writer.writerow(["varname", "value", "optid", "error"])
156        for name in PARAM_NAMES:
157            err = errors.get(name, 0) if errors else 0
158            writer.writerow([name, params[name], 0 if name in fix_set else 1, err])
159        writer.writerow(["VOC_est", params.get("VOC", 0), 0, 0])
160        writer.writerow(["ISC_est", params.get("ISC", 0), 0, 0])
161        if rss is not None: writer.writerow(["RSS_logy", rss, 0, 0])
162
163# --- 物理モデル・計算 ---
164
165def solve_root_poly(x, y, target_y=0, order=3):
166    """
167    多項式近似を用いて、与えられたデータセット(x, y)からyが特定の`target_y`になるxを推定します。
168
169    与えられたIVデータに対して、n次多項式で近似し、`y - target_y = 0`となる根を探します。
170    複数の実数根が見つかった場合、データの平均x値に最も近い根を選択します。
171    実数根が見つからない場合は、線形補間を使用します。
172
173    :param x: x軸のデータポイント（電圧）。
174    :type x: np.ndarray
175    :param y: y軸のデータポイント（電流）。
176    :type y: np.ndarray
177    :param target_y: 求めるyの値。デフォルトは0（開回路電圧VOCの推定に相当）。
178    :type target_y: float
179    :param order: 多項式近似の次数。デフォルトは3。
180    :type order: int
181    :returns: `target_y` に対応する推定されたx値。
182    :rtype: float
183    """
184    if len(x) < order + 1: order = max(1, len(x) - 1)
185    z = np.polyfit(x, y - target_y, order)
186    roots = np.roots(z)
187    real_roots = roots[np.isreal(roots)].real
188    if len(real_roots) == 0: return np.interp(target_y, y, x)
189    return real_roots[np.argmin(np.abs(real_roots - np.mean(x)))]
190
191def model(V, I0, ndiode, IPV, Rs, Rsh, temperature=300.0):
192    """
193    一ダイオードモデル（Single-Diode Model, SDM）に基づいてIV特性の電流を計算します。
194
195    このモデルは以下の数式で表されます：
196    I = IPV - I0 * (exp((V + I*Rs) / (ndiode * k * T / q)) - 1) - (V + I*Rs) / Rsh
197    ここで、Iは全電流、IPVは光電流、I0は飽和電流、ndiodeはダイオード理想因子、
198    Rsは直列抵抗、Rshはシャント抵抗、kはボルツマン定数、Tは温度、qは素電荷です。
199    この式はIについて解析的に解けないため、`scipy.optimize.brentq` を用いて数値的に解を求めます。
200    
201    :param V: 電圧のNumpy配列または単一値。
202    :type V: np.ndarray or float
203    :param I0: 飽和電流 (A)。
204    :type I0: float
205    :param ndiode: ダイオード理想因子。
206    :type ndiode: float
207    :param IPV: 光電流 (A)。
208    :type IPV: float
209    :param Rs: 直列抵抗 (Ω)。
210    :type Rs: float
211    :param Rsh: シャント抵抗 (Ω)。
212    :type Rsh: float
213    :param temperature: セル温度 (K)。デフォルトは300.0。
214    :type temperature: float
215    :returns: 入力電圧Vに対応する計算された電流 (A) のNumpy配列。
216    :rtype: np.ndarray
217    """
218    V = np.atleast_1d(V)
219    vt = ndiode * KB * temperature / E_CHARGE
220    def resid(Idark, v_val):
221        vd = v_val - Rs * Idark
222        arg = np.clip(vd / vt, -700, 700) # オーバーフロー防止
223        return I0 * (np.exp(arg) - 1.0) + vd / Rsh - Idark
224    results = []
225    for v_val in V:
226        try:
227            # 探索範囲を適切に設定
228            limit = abs(v_val / max(Rsh, EPS_R)) + abs(I0) * 1e5 + abs(IPV) + 1.0
229            id_val = brentq(resid, -limit, limit, args=(v_val,))
230            results.append(id_val - IPV)
231        except (ValueError, RuntimeError): # brentqが根を見つけられない場合
232            # フォールバックとして、Rsが非常に大きい場合の近似を返す
233            # または、単純な抵抗モデルに近い値を返す
234            results.append(v_val / max(Rsh, EPS_R) - IPV)
235    return np.array(results)
236
237def objective(p_free, V, I_meas, temp, base_params, fix_set):
238    """
239    最適化アルゴリズムのための目的関数を計算します。
240
241    この関数は、自由パラメータ`p_free`を更新し、それらの値を用いてモデルを呼び出し、
242    測定された電流と計算された電流の対数絶対値の差の二乗和（残差平方和）を返します。
243    `ndiode`を除くパラメータは対数スケールで最適化されます。
244    `EPS_I` は、対数計算でのゼロ除算や数値不安定性を防ぐために使用されます。
245
246    :param p_free: 最適化中の自由パラメータの配列。
247    :type p_free: np.ndarray
248    :param V: 測定された電圧データ。
249    :type V: np.ndarray
250    :param I_meas: 測定された電流データ。
251    :type I_meas: np.ndarray
252    :param temp: セル温度 (K)。
253    :type temp: float
254    :param base_params: 固定パラメータを含む、全てのパラメータの初期値または現在の値の辞書。
255    :type base_params: dict
256    :param fix_set: 固定されたパラメータ名のセット。
257    :type fix_set: set
258    :returns: 測定値とモデル予測値の対数絶対値の差の二乗和。
259    :rtype: float
260    """
261    curr = dict(base_params); idx = 0
262    for name in PARAM_NAMES:
263        if name not in fix_set:
264            val = p_free[idx]
265            curr[name] = 10**val if name != "ndiode" else val # ndiode以外はlog10スケールで最適化
266            idx += 1
267    I_calc = model(V, **{k: curr[k] for k in PARAM_NAMES}, temperature=temp)
268    # 微小値EPS_Iを加えてログ計算でのゼロ除算を防ぐ
269    return np.sum((np.log10(np.abs(I_meas) + EPS_I) - np.log10(np.abs(I_calc) + EPS_I))**2)
270
271# --- 誤差推定関連 ---
272
273def get_jacobian(p_free, V, temp, base_params, fix_set):
274    """
275    モデルの対数電流に対するヤコビ行列を数値的に計算します。
276
277    ヤコビ行列は、各自由パラメータに関するモデル出力（対数電流）の偏微分を近似したものです。
278    有限差分法を用いて計算され、誤差推定に使用されます。
279    `EPS_I` は、対数計算でのゼロ除算や数値不安定性を防ぐために使用されます。
280
281    :param p_free: 現在の自由パラメータの配列。
282    :type p_free: np.ndarray
283    :param V: 測定された電圧データ。
284    :type V: np.ndarray
285    :param temp: セル温度 (K)。
286    :type temp: float
287    :param base_params: 固定パラメータを含む、全てのパラメータの初期値または現在の値の辞書。
288    :type base_params: dict
289    :param fix_set: 固定されたパラメータ名のセット。
290    :type fix_set: set
291    :returns: モデルの対数電流に対するヤコビ行列。
292    :rtype: np.ndarray
293    """
294    eps = 1e-5 # 微小な摂動値
295    n_params, n_points = len(p_free), len(V)
296    J = np.zeros((n_points, n_params))
297
298    def get_log_model(p):
299        """与えられたパラメータでモデル電流を計算し、その対数絶対値を返すヘルパー関数"""
300        curr = dict(base_params); idx = 0
301        for name in PARAM_NAMES:
302            if name not in fix_set:
303                val = p[idx]
304                curr[name] = 10**val if name != "ndiode" else val
305                idx += 1
306        I_c = model(V, **{k: curr[k] for k in PARAM_NAMES}, temperature=temp)
307        return np.log10(np.abs(I_c) + EPS_I)
308    
309    f0 = get_log_model(p_free)
310    for i in range(n_params):
311        p_eps = np.copy(p_free); p_eps[i] += eps # 各パラメータを微小に摂動
312        J[:, i] = (get_log_model(p_eps) - f0) / eps # 有限差分で偏微分を計算
313    return J
314
315def estimate_errors(res, V, I_meas, temp, base_params, fix_set):
316    """
317    非線形最小二乗法で得られた結果に基づき、パラメータの誤差とモデル曲線の標準偏差を推定します。
318
319    フィッティング結果（`res`）とヤコビ行列を用いて、線形近似に基づいた共分散行列を計算します。
320    この共分散行列から、各パラメータの標準誤差と、モデル曲線自体の信頼性区間を推定します。
321    `np.linalg.pinv` を使用することで、数値的に不安定な場合でも計算を継続します。
322    パラメータの誤差は、`ndiode`は線形スケール、それ以外は元のスケール（対数スケールから変換）で報告されます。
323
324    :param res: `scipy.optimize.minimize` の結果オブジェクト。最適化されたパラメータ (`res.x`) と目的関数値 (`res.fun`) を含む。
325    :type res: scipy.optimize.OptimizeResult
326    :param V: 測定された電圧データ。
327    :type V: np.ndarray
328    :param I_meas: 測定された電流データ。
329    :type I_meas: np.ndarray
330    :param temp: セル温度 (K)。
331    :type temp: float
332    :param base_params: 固定パラメータを含む、全てのパラメータの初期値または現在の値の辞書。
333    :type base_params: dict
334    :param fix_set: 固定されたパラメータ名のセット。
335    :type fix_set: set
336    :returns:
337        - errors_dict (dict): 各自由パラメータの推定誤差の辞書。
338        - sigma_log_I (np.ndarray): モデル曲線の対数スケールでの標準偏差。
339    :rtype: tuple[dict, np.ndarray]
340    """
341    p_opt = res.x
342    n, p = len(I_meas), len(p_opt)
343    if n <= p: # データポイント数がパラメータ数以下の場合、誤差推定は不可能
344        return {}, np.zeros(len(V))
345    
346    J = get_jacobian(p_opt, V, temp, base_params, fix_set)
347    sigma2_hat = res.fun / (n - p) # 残差の分散を推定
348    
349    try:
350        # inv ではなく pinv (擬似逆行列) を使うことで、数値的に不安定な場合でも計算を継続
351        # J.T @ J は Hessian の近似（Gauss-Newton法における）
352        cov = sigma2_hat * np.linalg.pinv(J.T @ J) 
353        p_errors_log = np.sqrt(np.maximum(np.diag(cov), 0)) # パラメータの対数スケールでの標準誤差
354        sigma_log_I = np.sqrt(np.maximum(np.diag(J @ cov @ J.T), 0)) # モデル曲線の対数スケールでの標準偏差
355        
356        errors_dict, idx = {}, 0
357        for name in PARAM_NAMES:
358            if name not in fix_set:
359                if name == "ndiode": 
360                    errors_dict[name] = p_errors_log[idx] # ndiodeは線形スケール
361                else: 
362                    # log10スケールの誤差から元のスケールでの誤差に変換
363                    # d(10^x)/dx = 10^x * ln(10) を利用
364                    errors_dict[name] = (10**p_opt[idx]) * np.log(10) * p_errors_log[idx]
365                idx += 1
366        return errors_dict, sigma_log_I
367    except Exception as e:
368        print(f"Warning: Error estimation failed ({e})")
369        return {}, np.zeros(len(V))
370
371
372# --- 描画ロジック ---
373
374def plot_iv(V, I_meas, I_final, axes, sigma_log=None, label="Model", title=""):
375    """
376    測定IV特性とモデルIV特性を線形スケールおよび対数スケールでプロットします。
377
378    線形スケールのプロットと対数スケールのプロットの2つのサブプロット（axes）に、
379    測定データ、フィッティングされたモデル曲線、および（指定された場合）信頼性区間を描画します。
380    対数スケールプロットでは電流の絶対値が使用され、`EPS_I` を加えることでゼロ値のプロ処理を避けます。
381
382    :param V: 電圧データ。
383    :type V: np.ndarray
384    :param I_meas: 測定された電流データ。
385    :type I_meas: np.ndarray
386    :param I_final: モデルによって計算された最終電流データ。
387    :type I_final: np.ndarray
388    :param axes: プロットに使用するmatplotlibの2つのサブプロット軸 (ax_lin, ax_log)。
389    :type axes: tuple[matplotlib.axes.Axes, matplotlib.axes.Axes]
390    :param sigma_log: モデル曲線の対数スケールでの標準偏差。提供された場合、信頼性区間が描画されます。
391    :type sigma_log: np.ndarray or None
392    :param label: モデル曲線の凡例ラベル。デフォルトは"Model"。
393    :type label: str
394    :param title: プロットのタイトル。デフォルトは空文字列。
395    :type title: str
396    :returns: なし
397    :rtype: None
398    """
399    ax_lin, ax_log = axes
400    ax_lin.plot(V, I_meas, "o", color="blue", markersize=4, alpha=0.3, label="Data")
401    ax_lin.plot(V, I_final, "-", color="orange", lw=2, label=label)
402    
403    # 信頼性区間の塗りつぶし (水色領域)
404    if sigma_log is not None:
405        upper = 10**(np.log10(np.abs(I_final) + EPS_I) + sigma_log)
406        lower = 10**(np.log10(np.abs(I_final) + EPS_I) - sigma_log)
407        # np.sign(I_final) を乗算することで、電流の符号を維持しつつ信頼性区間をプロット
408        ax_lin.fill_between(V, np.sign(I_final)*lower, np.sign(I_final)*upper, color='lightblue', alpha=0.5)
409        # 対数プロットは絶対値でプロットされるため、符号は不要
410        ax_log.fill_between(V, lower, upper, color='lightblue', alpha=0.5)
411        
412    ax_lin.axhline(0, color='k', lw=0.5); ax_lin.axvline(0, color='k', lw=0.5)
413    ax_lin.set_title(title); ax_lin.grid(True); ax_lin.legend()
414    
415    # 対数プロットは絶対値で描画されるため、EPS_Iを加えてゼロを防ぐ
416    ax_log.semilogy(V, np.abs(I_meas) + EPS_I, "o", color="blue", markersize=4, alpha=0.3)
417    ax_log.semilogy(V, np.abs(I_final) + EPS_I, "-", color="orange", lw=2)
418    ax_log.grid(True, which="both")
419
420# --- モード別実行 ---
421
422def run_fit_mode(V, I_meas, init_params, args, csv_path, fix_set):
423    """
424    フィッティングモードを実行し、IV特性データに一ダイオードモデルを適合させます。
425
426    初期パラメータから最適化を開始し、`scipy.optimize.minimize` を用いて目的関数を最小化します。
427    フィッティングの進行状況は、指定された間隔でプロットとコンソール出力で表示されます。
428    最終的なフィッティングパラメータ、その誤差、および推定されたVOCとISCをCSVファイルに保存します。
429    `ndiode`を除くパラメータは対数スケールに変換されて最適化されます。
430
431    :param V: 測定された電圧データ。
432    :type V: np.ndarray
433    :param I_meas: 測定された電流データ。
434    :type I_meas: np.ndarray
435    :param init_params: フィッティングの初期パラメータの辞書。
436    :type init_params: dict
437    :param args: コマンドライン引数を格納したargparse.Namespaceオブジェクト。
438    :type args: argparse.Namespace
439    :param csv_path: 結果のパラメータを保存するCSVファイルのパス。
440    :type csv_path: str
441    :param fix_set: 固定するパラメータ名のセット。
442    :type fix_set: set
443    :returns:
444        - fit_params (dict): フィッティング後の最適化されたパラメータの辞書。
445        - sigma_log_curve (np.ndarray): 最適化されたモデル曲線の対数スケールでの標準偏差。
446    :rtype: tuple[dict, np.ndarray]
447    """
448    # ndiode以外はlog10スケールで最適化するため変換
449    p0_free = [init_params[p] if p == "ndiode" else np.log10(max(init_params[p], EPS_R if "R" in p else EPS_I)) 
450               for p in PARAM_NAMES if p not in fix_set]
451    
452    plt.ion(); fig, axes = plt.subplots(1, 2, figsize=(11, 5)); iters = [0]
453    
454    def callback(xk):
455        """最適化中の進捗をプロット・表示するためのコールバック関数"""
456        iters[0] += 1
457        if iters[0] % args.ninterval_plot == 0:
458            curr = dict(init_params); idx = 0
459            for name in PARAM_NAMES:
460                if name not in fix_set:
461                    curr[name] = 10**xk[idx] if name != "ndiode" else xk[idx]; idx += 1
462            I_c = model(V, **{k: curr[k] for k in PARAM_NAMES}, temperature=args.temperature)
463            rss = objective(xk, V, I_meas, args.temperature, init_params, fix_set)
464            for ax in axes: ax.cla() # 既存のプロットをクリア
465            plot_iv(V, I_meas, I_c, axes, label=f"Iter {iters[0]}", title=f"RSS: {rss:.2e}")
466            plt.pause(0.01) # 短時間一時停止してGUIを更新
467            if iters[0] % args.ninterval_print == 0:
468                p_str = " | ".join([f"{n}: {curr[n]:.2e}" + ("(F)" if n in fix_set else "") for n in PARAM_NAMES])
469                print(f"Iter {iters[0]:4d} | RSS: {rss:.4e} | {p_str}")
470
471    res = minimize(objective, p0_free, args=(V, I_meas, args.temperature, init_params, fix_set), 
472                   method=args.method, callback=callback)
473    
474    plt.ioff(); plt.close(fig) # インタラクティブモードを終了し、フィッティング中のプロットウィンドウを閉じる
475    
476    fit_params = dict(init_params); idx = 0
477    for name in PARAM_NAMES:
478        if name not in fix_set:
479            val = res.x[idx]; fit_params[name] = 10**val if name != "ndiode" else val; idx += 1
480    
481    errors, sigma_log_curve = estimate_errors(res, V, I_meas, args.temperature, init_params, fix_set)
482    
483    # 推定VOCとISCを計算
484    # ISCはV=0での電流値
485    fit_params["ISC"] = model(0, **{k: fit_params[k] for k in PARAM_NAMES}, temperature=args.temperature)[0]
486    # VOCはI=0となる電圧値（多項式根探索で推定）
487    v_f = np.linspace(min(V.min(), 0), max(V.max(), 1.2), 2000) # VOC推定のための広範囲電圧データ
488    fit_params["VOC"] = solve_root_poly(v_f, model(v_f, **{k: fit_params[k] for k in PARAM_NAMES}, temperature=args.temperature))
489    
490    save_param_csv(csv_path, fit_params, fix_set, errors=errors, rss=res.fun)
491    
492    print(f"\nFinal Results (RSS: {res.fun:.4e})")
493    for p in PARAM_NAMES:
494        err_str = f" +/- {errors[p]:.2e}" if p in errors else " (Fixed)"
495        print(f"  {p:8s}: {fit_params[p]:.6e}{err_str}")
496    return fit_params, sigma_log_curve
497
498def get_initial_params(V, I_meas, args, csv_path):
499    """
500    フィッティングのための初期パラメータを生成または既存のCSVファイルから読み込みます。
501
502    `--mode init` が指定されている場合、IV曲線の特徴（V=0近傍の電流からIPV、Vが負の領域の傾きからRsh）
503    から初期パラメータを推定します。
504    それ以外の場合、既存のパラメータCSVファイルから読み込みます。
505    CSVファイルが存在しない、または空の場合は、自動的に初期パラメータが推定されます。
506    コマンドライン引数で明示的に指定されたパラメータは、常に読み込み値や推定値を上書きします。
507
508    :param V: 測定された電圧データ。
509    :type V: np.ndarray
510    :param I_meas: 測定された電流データ。
511    :type I_meas: np.ndarray
512    :param args: コマンドライン引数を格納したargparse.Namespaceオブジェクト。
513    :type args: argparse.Namespace
514    :param csv_path: 初期パラメータを読み書きするCSVファイルのパス。
515    :type csv_path: str
516    :returns:
517        - params (dict): 初期パラメータ名とその値の辞書。
518        - csv_fix (set): CSVファイルで固定として指定されたパラメータ名のセット。
519    :rtype: tuple[dict, set]
520    """
521    if args.mode == "init":
522        # 初期パラメータの推定ロジック
523        # V=0近傍の電流からIPVを推定
524        uV, idxs = np.unique(V, return_index=True); uI = I_meas[idxs] # 重複電圧を除去
525        idx_v0 = np.argmin(np.abs(uV)) # V=0に最も近いインデックス
526        # V=0近傍のIV曲線から、直線近似/多項式近似でISCを推定（IPVに相当）
527        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)
528        # Vが負の領域の傾きからRshを推定
529        rsh_est = 1.0 / max(np.gradient(uI, uV)[uV < 0][0], 1e-12) if np.any(uV < 0) else 1e6
530        
531        params = {"I0": 1e-10, "ndiode": 1.5, "IPV": ipv_est, "Rs": 10.0, "Rsh": rsh_est, "VOC": 0.5, "ISC": -ipv_est}
532        csv_fix = set()
533    else:
534        params, csv_fix = load_param_csv(csv_path)
535        if not params: # CSVファイルが存在しない、または空の場合、initモードで初期化
536            print(f"Warning: Parameter CSV '{csv_path}' not found or empty. Generating initial parameters.")
537            temp_args = argparse.Namespace(**{**vars(args), 'mode': 'init'}) # 一時的にinitモードに設定
538            params, _ = get_initial_params(V, I_meas, temp_args, csv_path)
539    
540    # コマンドライン引数で指定されたパラメータで上書き
541    for p in PARAM_NAMES:
542        val = getattr(args, p); 
543        if val is not None: params[p] = val
544    return params, csv_fix
545
546# --- メイン ---
547
548def main():
549    """
550    スクリプトの主要な実行フローを制御します。
551
552    コマンドライン引数を解析し、指定されたモード（init, fit, sim）に基づいて、
553    データの読み込み、初期パラメータの取得、フィッティング実行、シミュレーション、
554    そして結果の保存とプロットを行います。
555    エラー発生時にはトレースバックを出力し、プログラムを終了します。
556    """
557    parser = argparse.ArgumentParser(description="一ダイオードモデルによる太陽電池IV特性フィッティング")
558    parser.add_argument("--mode", choices=["init", "fit", "sim"], default="fit",
559                        help="実行モード: init (初期パラメータ生成), fit (フィッティング), sim (シミュレーション)")
560    parser.add_argument("--method", default="Nelder-Mead",
561                        help="scipy.optimize.minimizeで使う最適化アルゴリズム (例: Nelder-Mead, L-BFGS-B)")
562    parser.add_argument("--infile", default="input.csv",
563                        help="入力IV特性CSVファイルのパス")
564    parser.add_argument("--temperature", type=float, default=300.0,
565                        help="計算に使用するセル温度 (K)")
566    parser.add_argument("--xmin", type=float,
567                        help="プロットおよびフィッティングに使用する電圧範囲の最小値 (V)")
568    parser.add_argument("--xmax", type=float,
569                        help="プロットおよびフィッティングに使用する電圧範囲の最大値 (V)")
570    for p in PARAM_NAMES: # 各モデルパラメータに対する引数を動的に追加
571        parser.add_argument(f"--{p}", type=float, help=f"パラメータ {p} の値を指定 (既存設定を上書き)")
572    parser.add_argument("--fix", nargs="*", default=[],
573                        help="フィッティング中に固定するパラメータの名前 (例: I0 ndiode)")
574    parser.add_argument("--nlsq_points", type=int, default=7, # 現在未使用だが、将来の拡張用
575                        help="非線形最小二乗近似のためのデータポイント数") 
576    parser.add_argument("--ninterval_print", type=int, default=10,
577                        help="フィッティング中のコンソール出力間隔 (イテレーション数)")
578    parser.add_argument("--ninterval_plot", type=int, default=10,
579                        help="フィッティング中のプロット更新間隔 (イテレーション数)")
580    args = parser.parse_args()
581
582    V, I_meas, inf = read_data(args.infile, xmin=args.xmin, xmax=args.xmax)
583    csv_path = Path(args.infile).stem + "-parameters.csv"
584    
585    params, csv_fix = get_initial_params(V, I_meas, args, csv_path)
586    fix_set = set(args.fix) | csv_fix # コマンドライン引数とCSV設定の両方で固定されたパラメータを結合
587    sigma_log = None
588
589    if args.mode == "init":
590        save_param_csv(csv_path, params, fix_set); print("\n--- Initial Parameters Saved ---")
591    elif args.mode == "fit":
592        params, sigma_log = run_fit_mode(V, I_meas, params, args, csv_path, fix_set)
593    elif args.mode == "sim":
594        print("\n--- Simulation Mode ---")
595        for p in PARAM_NAMES: print(f"  {p:8s}: {params[p]:.6e}")
596
597    # 最終結果のプロット
598    I_final = model(V, **{k: params[k] for k in PARAM_NAMES}, temperature=args.temperature)
599    fig, axes = plt.subplots(1, 2, figsize=(11, 5))
600    plot_iv(V, I_meas, I_final, axes, sigma_log=sigma_log, label=f"Final ({args.mode})", title=inf["FileName"])
601    plt.tight_layout(); plt.show()
602
603if __name__ == "__main__":
604    try: main()
605    except Exception: traceback.print_exc(); sys.exit(1)