#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
拡張一ダイオードモデル（SDM）IVフィッティングツール
このモジュールは、太陽電池や他の半導体デバイスの電流-電圧（IV）特性を、
拡張一ダイオードモデル（SDM）を用いてフィッティング、初期値推定、またはシミュレーションするためのツールを提供します。
複数の輸送メカニズム（ダイオード、トンネル誘起フォワードエミッション (TFE)、
非理想フォワード電流 (FN)、空間電荷制限電流 (SCLC)）を組み合わせて利用できます。
機能:
- フォワードバイアスとリバースバイアスのモデルを独立に選択可能
- ダイオード、TFE、FN、SCLCの各メカニズムを任意に組み合わせ可能
- 数値解法は接合電圧 Vd を未知数とする根探しで安定化
- モデルに未使用のパラメータは自動的に固定
- forwardIV と reverseIV が同一の場合、再計算を回避して効率化
- 最適化途中のIVカーブをアニメーション表示（常に有効）
- 線形近似に基づくパラメータ誤差推定
- モデル電流の標準偏差から信頼性区間を水色領域で表示
mode=init の初期値推定:
- 入力IVを電圧昇順にソート
- I-Vを局所多項式で平滑化
- dI/dV から Rs, Rsh を推定
- ISC = I(V=0) を V=0 近傍の多項式で推定
- IPV = -ISC
- VOC = I(V)=0 の多項式根で推定
- Ish = V<0 の代表点 Vsh における電流
- I0 = ISC - Ish （シャント電流は差し引かない）
- TFE 初期値:
V<0 の十分負側で ln(|I|)-V を一次近似
E00 = 1/|slope|
A_tfe = |I| * exp(+V/E00) を代表点から推定
- ndiode は固定初期値 1.5
- FN / SCLC の初期値は結果に影響しない極小値
注意:
- mode=init の初期値は fit 用の頑健な初期値です。
- 物理的に厳密な意味づけは mode=fit で精密化してください。
関連リンク:
:doc:`pvfit_usage`
"""
import os
import sys
import argparse
import builtins
from pathlib import Path
import traceback
import csv
from openpyxl import Workbook
from openpyxl.styles import Font, PatternFill
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import brentq, minimize
KB = 1.380649e-23
E_CHARGE = 1.602176634e-19
PARAM_NAMES = [
"I0", "ndiode", "IPV", "Rs", "Rsh",
"A_tfe", "E00",
"A_fn", "B_fn",
"K_sclc",
]
LOG_PARAMS = {
"I0", "IPV", "Rs", "Rsh",
"A_tfe", "E00",
"A_fn", "B_fn",
"K_sclc",
}
EPS_I = 1.0e-15
EPS_R = 1.0e-12
DEFAULT_PARAMETERS = {
"I0": 1.0e-10,
"ndiode": 1.5,
"IPV": 1.0e-6,
"Rs": 1.0,
"Rsh": 1.0e9,
"A_tfe": 1.0e-30,
"E00": 1.0e-3,
"A_fn": 1.0e-30,
"B_fn": 10.0,
"K_sclc": 1.0e-30,
"VOC_est": 0.5,
"ISC_est": -1.0e-6,
"Vsh_est": -0.1,
"Ish_est": -2.0e-6,
"Vtfe_est": -0.2,
"Itfe_est": 1.0e-12,
}
fontsize = 16
[ドキュメント]
def initialize():
"""
コマンドライン引数パーサーを初期化し、引数を解析します。
各種実行モード (`analyze`, `init`, `fit`, `sim`)、最適化手法、
入力ファイル、温度、データクリッピング範囲、スキップ間隔、
フォワード/リバースIVモデル、ブレンド電圧幅、プロット/プリント間隔、
最大イテレーション数、許容誤差、固定パラメータなどの引数を定義します。
:returns: args (解析された引数を含むオブジェクト), parser (ArgumentParserオブジェクト)
:rtype: tuple[argparse.Namespace, argparse.ArgumentParser]
"""
parser = argparse.ArgumentParser(description="Extended SDM IV fitting tool with diode/TFE/FN/SCLC.")
parser.add_argument("--mode", choices=["analyze", "init", "fit", "sim"], default="fit",
help="Execution mode (init: initial parameter estimation, fit: fit parameters to data, sim: simulate IV curve)")
parser.add_argument("--method", default="Nelder-Mead",
help="scipy.optimize.minimize method (e.g., 'Nelder-Mead', 'Powell', 'BFGS', 'L-BFGS-B')")
parser.add_argument("--infile", help="Input IV CSV file path")
parser.add_argument("--temperature", type=float, default=300.0,
help="Temperature in K for the model")
parser.add_argument("--xmin", type=float, help="Minimum voltage (V) for data clipping, fitting, and plotting")
parser.add_argument("--xmax", type=float, help="Maximum voltage (V) for data clipping, fitting, and plotting")
parser.add_argument("--ndataskip", type=int, default=0, help="# of data to skip reading")
parser.add_argument("--forwardIV", default="diode",
help="Transport mechanisms for forward-bias side, e.g., 'diode' or 'diode+tfe+sclc'")
parser.add_argument("--reverseIV", default="diode+tfe",
help="Transport mechanisms for reverse-bias side, e.g., 'diode+tfe' or 'diode+fn'")
parser.add_argument("--dv_blend", type=float, default=0.05,
help="Sigmoid transition width (V) between reverse and forward models for junction voltage Vd")
parser.add_argument("--ninterval_print", type=int, default=10,
help="Console print interval during fitting iterations")
parser.add_argument("--ninterval_plot", type=int, default=10,
help="Animation update interval during fitting iterations")
parser.add_argument("--nmaxiter", type=int, default=1000,
help="Maximum optimizer iterations")
parser.add_argument("--tol", type=float, default=1.0e-7,
help="Optimizer tolerance (xatol and fatol)")
parser.add_argument("--fix", nargs="*", default=[],
help="Space-separated list of parameter names to fix during fitting (e.g., 'I0 ndiode')")
parser.add_argument("--outprefix", default="pvfit", help="Output file prefix")
for p in PARAM_NAMES:
parser.add_argument(f"--{p}", type=float, help=f"Override initial value for {p}")
args = parser.parse_args()
args.animate_fit = True # 最適化中のアニメーション表示は常に有効
return args, parser
# =============================================================================
# argument repeat / validation
# =============================================================================
[ドキュメント]
def print_args_and_derived(args):
"""
コマンドライン引数とその派生値をコンソールに出力します。
ユーザーが指定した入力値と、プログラム内部で決定される追加の設定（例：使用されるパラメータ、自動固定されるパラメータ）
を一覧表示し、設定の確認を容易にします。
:param args: 解析されたコマンドライン引数を含むオブジェクト。
:type args: argparse.Namespace
:returns: None
:rtype: None
"""
print("=== Arguments ===")
print(f"mode = {args.mode}")
print(f"infile = {args.infile}")
print(f"temperature = {args.temperature} K")
print(f"xmin = {args.xmin} V" if args.xmin is not None else "xmin = None")
print(f"xmax = {args.xmax} V" if args.xmax is not None else "xmax = None")
print(f"ndataskip = {args.ndataskip}")
print(f"forwardIV = {args.forwardIV}")
print(f"reverseIV = {args.reverseIV}")
print(f"dv_blend = {args.dv_blend} V")
print(f"ninterval_print = {args.ninterval_print}")
print(f"ninterval_plot = {args.ninterval_plot}")
print(f"nmaxiter = {args.nmaxiter}")
print(f"tol = {args.tol}")
print(f"fix = {' '.join(args.fix) if args.fix else '(none)'}")
print(f"outprefix = {args.outprefix}")
for name in PARAM_NAMES:
val = getattr(args, name, None)
print(f"{name:12s} = {val}" if val is not None else f"{name:12s} = None")
print()
print("=== Derived values ===")
used = used_params_from_modes(args.forwardIV, args.reverseIV)
auto_fixed = set(PARAM_NAMES) - used
print(f"used_params = {' '.join(sorted(used)) if used else '(none)'}")
print(f"auto_fixed = {' '.join(sorted(auto_fixed)) if auto_fixed else '(none)'}")
print()
[ドキュメント]
def validate_args(args):
"""
コマンドライン引数の有効性を検証します。
モードごとの必須引数、数値引数の範囲、パラメータ名の整合性などをチェックし、
無効な引数が存在する場合は `ValueError` を発生させます。
:param args: 解析されたコマンドライン引数を含むオブジェクト。
:type args: argparse.Namespace
:returns: None
:rtype: None
:raises ValueError: 無効な引数が見つかった場合。
"""
if args.mode in ("init", "fit", "analyze"):
if not args.infile:
raise ValueError(f"--infile is required for mode={args.mode}")
if args.ndataskip < 0:
raise ValueError(f"ndataskip must be >= 0, got {args.ndataskip}")
if args.xmin is not None and args.xmax is not None and args.xmin > args.xmax:
raise ValueError(f"xmin must be <= xmax, got xmin={args.xmin}, xmax={args.xmax}")
if args.temperature <= 0.0:
raise ValueError(f"temperature must be > 0 K, got {args.temperature}")
if args.dv_blend < 0.0:
raise ValueError(f"dv_blend must be >= 0 V, got {args.dv_blend}")
if args.ninterval_print <= 0:
raise ValueError(f"ninterval_print must be > 0, got {args.ninterval_print}")
if args.ninterval_plot <= 0:
raise ValueError(f"ninterval_plot must be > 0, got {args.ninterval_plot}")
if args.nmaxiter <= 0:
raise ValueError(f"nmaxiter must be > 0, got {args.nmaxiter}")
if args.tol <= 0.0:
raise ValueError(f"tol must be > 0, got {args.tol}")
unknown_fix = [x for x in args.fix if x not in PARAM_NAMES]
if unknown_fix:
raise ValueError(f"Unknown parameter name(s) in --fix: {unknown_fix}")
for name in LOG_PARAMS:
val = getattr(args, name, None)
if val is not None and val <= 0.0:
raise ValueError(f"{name} must be > 0 because it is log-scaled, got {val}")
# 非logパラメータ側
if getattr(args, "ndiode", None) is not None and args.ndiode <= 0.0:
raise ValueError(f"ndiode must be > 0, got {args.ndiode}")
if getattr(args, "B_fn", None) is not None and args.B_fn <= 0.0:
raise ValueError(f"B_fn must be > 0, got {args.B_fn}")
# =============================================================================
# data I/O
# =============================================================================
[ドキュメント]
def read_data(infile, xmin=None, xmax=None, ndataskip=0):
"""
指定されたCSVファイルからIVデータ（電圧と電流）を読み込みます。
ファイルからメタデータを抽出し、電圧昇順にデータをソートし、
指定された電圧範囲でクリッピングし、データスキップ間隔を適用します。
:param infile: 入力IV CSVファイルのパス。
:type infile: str
:param xmin: データをクリップする最小電圧 (V)。Noneの場合、下限なし。
:type xmin: float or None
:param xmax: データをクリップする最大電圧 (V)。Noneの場合、上限なし。
:type xmax: float or None
:param ndataskip: 読み込み時にスキップするデータ点の数（n+1点ごとにデータを保持）。
:type ndataskip: int
:returns: V (読み込まれた電圧データのNumPy配列), I (読み込まれた電流データのNumPy配列), info (ファイル名と記録時間を含む辞書)
:rtype: tuple[numpy.ndarray, numpy.ndarray, dict]
:raises ValueError: 有効なデータが見つからない場合、またはクリッピング後にデータ点が残らない場合。
"""
data_started = False
raw_v = []
raw_i = []
info = {"FileName": infile, "RecordTime": "Unknown"}
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) >= 3 and cols[0] == "MetaData" and "RecordTime" in cols[1]:
info["RecordTime"] = cols[2]
if cols[0] == "DataName":
data_started = True
continue
if data_started and cols[0] == "DataValue":
if len(cols) < 3:
continue
v = float(cols[1])
i = float(cols[2])
if not np.isfinite(v) or not np.isfinite(i):
continue
raw_v.append(v)
raw_i.append(i)
if not raw_v:
raise ValueError(f"No valid data found in {infile}")
V = np.asarray(raw_v, dtype=float)
I = np.asarray(raw_i, dtype=float)
# まず V でソート
idx = np.argsort(V)
V = V[idx]
I = I[idx]
# 次に xmin/xmax でクリップ
mask = np.ones(len(V), dtype=bool)
if xmin is not None:
mask &= (V >= xmin)
if xmax is not None:
mask &= (V <= xmax)
V = V[mask]
I = I[mask]
if len(V) == 0:
raise ValueError("No data points remain after xmin/xmax clipping.")
# 最後に ndataskip を適用
if ndataskip > 0:
keep = np.zeros(len(V), dtype=bool)
keep[::(ndataskip + 1)] = True
V = V[keep]
I = I[keep]
print(f"File Name : {info['FileName']}")
print(f"Record Time : {info['RecordTime']}")
print(f"Points : {len(V)}")
print(f"V range : {V.min():.6g} to {V.max():.6g} V")
return V, I, info
[ドキュメント]
def load_param_csv(csv_path):
"""
パラメータ設定CSVファイルからパラメータ値と固定設定を読み込みます。
CSVファイルには 'varname', 'value', 'optid' の列が含まれると想定されます。
'optid' が '0' の場合、そのパラメータは固定されます。
:param csv_path: 読み込むCSVファイルのパス。
:type csv_path: str
:returns: params (読み込まれたパラメータ名と値の辞書), fix_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.get("varname", "").strip()
if not name:
continue
try:
val = float(row.get("value", ""))
except (TypeError, ValueError):
continue
params[name] = val
try:
if int(row.get("optid", 1)) == 0 and name in PARAM_NAMES:
fix_set.add(name)
except (TypeError, ValueError):
pass
return params, fix_set
[ドキュメント]
def save_param_csv(csv_path, params, fix_set, errors=None, rss=None):
"""
現在のパラメータ値、固定設定、推定誤差、およびRSSをCSVファイルに保存します。
`PARAM_NAMES` にリストされている全てのパラメータと、追加の診断情報（VOC_estなど）が保存されます。
:param csv_path: 保存するCSVファイルのパス。
:type csv_path: str
:param params: 保存するパラメータ名と値の辞書。
:type params: dict
:param fix_set: 固定されたパラメータ名のセット。
:type fix_set: set
:param errors: 各パラメータの推定誤差を含む辞書。Noneの場合、誤差は0.0として保存。
:type errors: dict or None
:param rss: 残差平方和 (Residual Sum of Squares)。Noneの場合、保存しない。
:type rss: float or None
:returns: None
: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:
writer.writerow([
name,
params[name],
0 if name in fix_set else 1,
errors.get(name, 0.0),
])
for extra in [
"VOC_est", "ISC_est", "Vsh_est", "Ish_est",
"Vtfe_est", "Itfe_est",
]:
if extra in params:
writer.writerow([extra, params[extra], 0, 0.0])
if rss is not None:
writer.writerow(["RSS_logy", rss, 0, 0.0])
[ドキュメント]
def save_iv_to_excel(outfile_xlsx, V, I_input, I_init, I_fit, info,
params_init=None, params_final=None, rss=None, errors=None):
"""
IVデータ、初期パラメータ、最終パラメータ、およびサマリー情報をExcelファイルに保存します。
測定データ、初期モデル電流、フィッティング後のモデル電流が線形および対数絶対値スケールで保存されます。
パラメータシートには初期値と最終値、および推定誤差が含まれます。
:param outfile_xlsx: 出力するExcelファイルのパス。
:type outfile_xlsx: str
:param V: 電圧データのNumPy配列。
:type V: numpy.ndarray
:param I_input: 測定された電流データのNumPy配列。Noneの場合、保存しない。
:type I_input: numpy.ndarray or None
:param I_init: 初期モデルによって計算された電流データのNumPy配列。Noneの場合、保存しない。
:type I_init: numpy.ndarray or None
:param I_fit: フィッティング後のモデルによって計算された電流データのNumPy配列。Noneの場合、保存しない。
:type I_fit: numpy.ndarray or None
:param info: ファイル名や記録時間などの情報を含む辞書。
:type info: dict
:param params_init: 初期パラメータ値を含む辞書。Noneの場合、保存しない。
:type params_init: dict or None
:param params_final: フィッティング後の最終パラメータ値を含む辞書。Noneの場合、保存しない。
:type params_final: dict or None
:param rss: 残差平方和 (Residual Sum of Squares)。Noneの場合、保存しない。
:type rss: float or None
:param errors: 各パラメータの推定誤差を含む辞書。Noneの場合、保存しない。
:type errors: dict or None
:returns: None
:rtype: None
"""
print(f"[I/O] Save IV Excel: {outfile_xlsx}")
wb = Workbook()
header_fill = PatternFill("solid", fgColor="1F4E78")
header_font = Font(color="FFFFFF", bold=True)
def log_abs_current(x):
if x is None:
return None
return np.log10(np.abs(x) + EPS_I)
# ------------------------------------------------------------
# sheet 1: iv_data
# ------------------------------------------------------------
ws = wb.active
ws.title = "iv_data"
headers = [
"V",
"I_input", "I_init", "I_fit",
"abs_I_input", "abs_I_init", "abs_I_fit",
"log10_abs_I_input", "log10_abs_I_init", "log10_abs_I_fit",
]
ws.append(headers)
for c in ws[1]:
c.fill = header_fill
c.font = header_font
n = len(V)
for i in range(n):
ii = None if I_input is None else float(I_input[i])
i0 = None if I_init is None else float(I_init[i])
ifit = None if I_fit is None else float(I_fit[i])
ws.append([
float(V[i]),
ii, i0, ifit,
None if ii is None else abs(ii),
None if i0 is None else abs(i0),
None if ifit is None else abs(ifit),
None if ii is None else float(log_abs_current(ii)),
None if i0 is None else float(log_abs_current(i0)),
None if ifit is None else float(log_abs_current(ifit)),
])
for col in ["A","B","C","D","E","F","G","H","I","J"]:
ws.column_dimensions[col].width = 18
# ------------------------------------------------------------
# sheet 2: init_params
# ------------------------------------------------------------
if params_init is not None:
ws2 = wb.create_sheet("init_params")
ws2.append(["name", "value"])
for c in ws2[1]:
c.fill = header_fill
c.font = header_font
for name in PARAM_NAMES:
ws2.append([name, float(params_init.get(name, np.nan))])
for extra in ["VOC_est", "ISC_est", "Vsh_est", "Ish_est", "Vtfe_est", "Itfe_est"]:
if extra in params_init:
ws2.append([extra, float(params_init[extra])])
# ------------------------------------------------------------
# sheet 3: final_params
# ------------------------------------------------------------
if params_final is not None:
ws3 = wb.create_sheet("final_params")
ws3.append(["name", "value", "error"])
for c in ws3[1]:
c.fill = header_fill
c.font = header_font
if errors is None:
errors = {}
for name in PARAM_NAMES:
ws3.append([
name,
float(params_final.get(name, np.nan)),
float(errors.get(name, np.nan)),
])
# ------------------------------------------------------------
# sheet 4: summary
# ------------------------------------------------------------
ws4 = wb.create_sheet("summary")
ws4.append(["item", "value"])
for c in ws4[1]:
c.fill = header_fill
c.font = header_font
ws4.append(["FileName", info.get("FileName", "")])
ws4.append(["RecordTime", info.get("RecordTime", "")])
ws4.append(["Npoints", int(len(V))])
if rss is not None:
ws4.append(["RSS_logy", float(rss)])
wb.save(outfile_xlsx)
print(f"[I/O] Saved IV Excel: {outfile_xlsx}")
# =============================================================================
# helper math
# =============================================================================
[ドキュメント]
def solve_root_poly(x, y, target_y=0.0, order=3):
"""
データ点 (x, y) に多項式をフィッティングし、指定された `target_y` に対応する x の根を近似的に見つけます。
データ点数が少ない場合は、線形補間または多項式の次数を調整して対処します。
複数の実根がある場合は、入力 `x` の平均に最も近い根を選択します。
:param x: 独立変数データ点の配列。
:type x: numpy.ndarray
:param y: 従属変数データ点の配列。
:type y: numpy.ndarray
:param target_y: 根を見つけたい目標のy値。
:type target_y: float
:param order: 多項式フィッティングの次数。
:type order: int
:returns: `target_y` に対応する x の推定根。
:rtype: float
"""
if len(x) < 2:
return float(x[0]) if len(x) else 0.0
order = min(order, len(x) - 1)
if order < 1:
return float(np.interp(target_y, y, x))
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 float(np.interp(target_y, y, x))
return float(real_roots[np.argmin(np.abs(real_roots - np.mean(x)))])
[ドキュメント]
def local_poly_value(x, y, x0, npts=7, order=3):
"""
指定された点 `x0` の近傍のデータ点に局所多項式をフィッティングし、`x0` におけるy値を推定します。
`npts` で指定された数の最も近いデータ点を使用します。
:param x: 独立変数データ点の配列。
:type x: numpy.ndarray
:param y: 従属変数データ点の配列。
:type y: numpy.ndarray
:param x0: y値を推定したい目標のx値。
:type x0: float
:param npts: 局所フィッティングに使用するデータ点の数。
:type npts: int
:param order: 多項式フィッティングの次数。
:type order: int
:returns: `x0` における推定y値。
:rtype: float
"""
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
idx = np.argsort(np.abs(x - x0))[:max(2, npts)]
xs = x[idx]
ys = y[idx]
sidx = np.argsort(xs)
xs = xs[sidx]
ys = ys[sidx]
deg = min(order, len(xs) - 1)
if deg < 1: # データ点が1点のみの場合、その点のy値を返す
if len(xs) == 1:
return float(ys[0])
else: # データ点が2点の場合、線形補間
return float(np.interp(x0, xs, ys))
coef = np.polyfit(xs, ys, deg)
return float(np.polyval(coef, x0))
[ドキュメント]
def smooth_polyfit(y, window_points=5, poly_order=3):
"""
Savitzky-Golayフィルターに似た局所多項式フィッティングにより、配列 `y` を平滑化します。
各データ点に対して、その点と近傍の `window_points` 数（奇数）のデータ点に
`poly_order` 次の多項式をフィッティングし、中央の点の値を推定します。
:param y: 平滑化するデータ点の配列。
:type y: numpy.ndarray
:param window_points: 各局所フィッティングに使用するウィンドウの点の数（奇数であるべき）。
:type window_points: int
:param poly_order: 各局所フィッティングで使用する多項式の次数。
:type poly_order: int
:returns: 平滑化されたデータ点の配列。
:rtype: numpy.ndarray
"""
y = np.asarray(y, dtype=float)
n = len(y)
if n < 3:
return y.copy()
if window_points < 3:
window_points = 3
if window_points % 2 == 0:
window_points += 1
if window_points > n:
window_points = n if (n % 2 == 1) else n - 1
if window_points < 3: # ここに到達したら、平滑化不能なのでコピーを返す
return y.copy()
if poly_order >= window_points:
poly_order = window_points - 1
if poly_order < 1:
poly_order = 1
half = window_points // 2
ys = np.empty(n, dtype=float)
for i in range(n):
i0 = max(0, i - half)
i1 = min(n, i + half + 1)
while (i1 - i0) < window_points:
if i0 > 0:
i0 -= 1
elif i1 < n:
i1 += 1
else:
break
idx = np.arange(i0, i1, dtype=float)
yy = y[i0:i1]
xloc = idx - i # 中心を0とする相対座標
deg = min(poly_order, len(yy) - 1)
if deg < 1:
ys[i] = y[i]
continue
coeff = np.polyfit(xloc, yy, deg)
ys[i] = np.polyval(coeff, 0.0)
return ys
[ドキュメント]
def sigmoid_blend(x, dV):
"""
シグモイド関数を用いて、2つのモデル間のブレンド係数を計算します。
`x` が0.0の場合に0.5、`x` が正で十分に大きい場合に1.0、
`x` が負で十分に小さい場合に0.0に近づく値を返します。
`dV` はブレンドの遷移幅を決定します。`dV` が非常に小さい場合、ステップ関数に近づきます。
:param x: ブレンドの基準となる電圧または汎用的な値。
:type x: numpy.ndarray or float
:param dV: シグモイド遷移の幅を制御するパラメータ。
:type dV: float
:returns: `x` に応じたブレンド係数。0から1の間の値を取ります。
:rtype: numpy.ndarray or float
"""
x = np.asarray(x, dtype=float)
if dV < 1.0e-6:
return (x >= 0.0).astype(float)
z = np.clip(-x / dV, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
return 1.0 / (1.0 + np.exp(z))
# =============================================================================
# current components
# =============================================================================
[ドキュメント]
def j_diode(Vd, I0, ndiode, temperature):
"""
理想ダイオードモデル（Shockley方程式）の電流成分を計算します。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param I0: 逆飽和電流。
:type I0: float
:param ndiode: ダイオード理想因子（エミッタ係数）。
:type ndiode: float
:param temperature: 温度 (K)。
:type temperature: float
:returns: ダイオード電流。
:rtype: numpy.ndarray or float
"""
vt = ndiode * KB * temperature / E_CHARGE
arg = np.clip(Vd / vt, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
return I0 * (np.exp(arg) - 1.0)
[ドキュメント]
def j_tfe_forward(Vd, A_tfe, E00):
"""
熱アシストトンネルフィールドエミッション (TFE) 電流のフォワードバイアス成分を計算します。
TFEは通常、低温や欠陥の多い材料で支配的になるメカニズムです。
この関数は、Vd > 0 の領域でTFEによる電流を計算します。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param A_tfe: TFE電流の振幅係数。
:type A_tfe: float
:param E00: TFEの特性エネルギーまたは係数。
:type E00: float
:returns: TFEフォワード電流。
:rtype: numpy.ndarray or float
"""
e00 = max(E00, 1.0e-12) # ゼロ除算防止
arg = np.clip(Vd / e00, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
return A_tfe * (np.exp(arg) - 1.0)
[ドキュメント]
def j_tfe_reverse(Vd, A_tfe, E00):
"""
熱アシストトンネルフィールドエミッション (TFE) 電流のリバースバイアス成分を計算します。
この関数は、Vd < 0 の領域でTFEによる電流を計算します。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param A_tfe: TFE電流の振幅係数。
:type A_tfe: float
:param E00: TFEの特性エネルギーまたは係数。
:type E00: float
:returns: TFEリバース電流。
:rtype: numpy.ndarray or float
"""
e00 = max(E00, 1.0e-12) # ゼロ除算防止
arg = np.clip(-Vd / e00, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
return -A_tfe * (np.exp(arg) - 1.0)
[ドキュメント]
def j_fn_forward(Vd, A_fn, B_fn):
"""
Fowler-Nordheimトンネル電流 (FN) のフォワードバイアス成分を計算します。
このメカニズムは、高注入レベルや特定の欠陥メカニズムによって引き起こされる電流をモデル化するために使用できます。
Vd > 0 の領域で定義されます。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param A_fn: FN電流の振幅係数。
:type A_fn: float
:param B_fn: FN電流の指数関数的電圧依存性を制御する係数。
:type B_fn: float
:returns: FNフォワード電流。
:rtype: numpy.ndarray or float
"""
if np.ndim(Vd) == 0:
if Vd <= 1.0e-12:
return 0.0
arg = np.clip(-B_fn / Vd, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
return A_fn * (Vd ** 2) * np.exp(arg)
arr = np.asarray(Vd, dtype=float)
out = np.zeros_like(arr)
mask = arr > 1.0e-12
arg = np.clip(-B_fn / arr[mask], -700.0, 700.0) # オーバーフロー防止のためのクリッピング
out[mask] = A_fn * (arr[mask] ** 2) * np.exp(arg)
return out
[ドキュメント]
def j_fn_reverse(Vd, A_fn, B_fn):
"""
Fowler-Nordheimトンネル電流 (FN) のリバースバイアス成分を計算します。
このメカニズムは、高注入レベルや特定の欠陥メカニズムによって引き起こされる電流をモデル化するために使用できます。
Vd < 0 の領域で定義されます。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param A_fn: FN電流の振幅係数。
:type A_fn: float
:param B_fn: FN電流の指数関数的電圧依存性を制御する係数。
:type B_fn: float
:returns: FNリバース電流。
:rtype: numpy.ndarray or float
"""
if np.ndim(Vd) == 0:
if Vd >= -1.0e-12:
return 0.0
vabs = -Vd
arg = np.clip(-B_fn / vabs, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
return -A_fn * (vabs ** 2) * np.exp(arg)
arr = np.asarray(Vd, dtype=float)
out = np.zeros_like(arr)
mask = arr < -1.0e-12
vabs = -arr[mask]
arg = np.clip(-B_fn / vabs, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
out[mask] = -A_fn * (vabs ** 2) * np.exp(arg)
return out
[ドキュメント]
def j_sclc_trap_transition(Vd, K_sclc, m_sclc, Vtfl, dV=1e-3):
"""
トラップ制限空間電荷制限電流（SCLC）とトラップフリーSCLCの間の遷移をモデル化する電流成分を計算します。
このモデルは、電圧に応じてトラップが満たされることで、電流が異なる電圧依存性を示す領域を表現します。
具体的には、トラップ制限領域では `Vd^m_sclc` に、トラップフリー領域では `Vd^2` に従います。
シグモイド関数を用いて両者間の滑らかな遷移を実現します。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param K_sclc: SCLCの振幅係数（トラップ制限領域の係数）。
:type K_sclc: float
:param m_sclc: トラップ制限領域でのSCLCの電圧依存性指数。
:type m_sclc: float
:param Vtfl: トラップが満たされる閾値電圧。
:type Vtfl: float
:param dV: 遷移の幅を制御する電圧（シグモイド関数の傾き）。
:type dV: float
:returns: トラップ遷移を考慮したSCLC電流。
:rtype: numpy.ndarray or float
"""
vabs = np.abs(Vd)
K_tf = K_sclc * (max(Vtfl, 1.0e-30) ** (m_sclc - 2.0))
if dV < 1.0e-12:
w = (vabs >= Vtfl).astype(float) if np.ndim(vabs) else float(vabs >= Vtfl)
else:
z = np.clip(-(vabs - Vtfl) / dV, -700.0, 700.0)
w = 1.0 / (1.0 + np.exp(z))
i_low = K_sclc * (vabs ** m_sclc) # trap-limited
i_high = K_tf * (vabs ** 2.0) # trap-free
imag = (1.0 - w) * i_low + w * i_high
return np.sign(Vd) * imag
[ドキュメント]
def j_sclc(Vd, K_sclc):
"""
空間電荷制限電流 (SCLC) を計算します。
SCLCは、高抵抗材料におけるキャリア注入が空間電荷によって制限されるときに発生します。
この電流は電圧の2乗に比例し、電圧の符号によって方向が決まります。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param K_sclc: SCLCの振幅係数。
:type K_sclc: float
:returns: SCLC電流。
:rtype: numpy.ndarray or float
"""
return np.sign(Vd) * K_sclc * (np.abs(Vd) ** 2)
[ドキュメント]
def parse_mech_mode(mode_str):
"""
輸送メカニズムのモード文字列を解析し、使用するメカニズムのセットを返します。
モード文字列はカンマ (`,`) またはプラス (`+`) で区切られたメカニズム名（例: "diode+tfe"）を
含みます。無効なメカニズム名が指定された場合はValueErrorを発生させます。
:param mode_str: 解析するモード文字列（例: "diode", "diode+tfe", "none"）。
:type mode_str: str or None
:returns: 使用するメカニズム名のセット。
:rtype: set
:raises ValueError: 未知の輸送メカニズム名が指定された場合。
"""
if mode_str is None:
return set()
txt = mode_str.lower().replace(" ", "")
if txt in ("", "none", "off"):
return set()
parts = [p for p in txt.replace(",", "+").split("+") if p]
norm = set()
for p in parts:
if p in ("diode", "tfe", "fn", "sclc"):
norm.add(p)
# elif p == "fn": # エイリアス: 既に上記で処理されているため、この行は冗長
# norm.add("fn")
else:
raise ValueError(f"Unknown transport mechanism: {p}")
return norm
[ドキュメント]
def used_params_from_modes(forwardIV, reverseIV):
"""
指定されたフォワードおよびリバースIVモデルで使用されるパラメータ名を決定します。
共通のパラメータ (`Rs`, `Rsh`, `IPV`) に加えて、
有効な輸送メカニズム (`diode`, `tfe`, `fn`, `sclc`) に関連するパラメータを特定します。
:param forwardIV: フォワードバイアスモデルのメカニズムを示す文字列（例: "diode", "diode+tfe"）。
:type forwardIV: str
:param reverseIV: リバースバイアスモデルのメカニズムを示す文字列（例: "diode+tfe", "diode+fn"）。
:type reverseIV: str
:returns: 使用されるパラメータ名のセット。
:rtype: set
"""
used = {"Rs", "Rsh", "IPV"}
all_modes = parse_mech_mode(forwardIV) | parse_mech_mode(reverseIV)
if "diode" in all_modes:
used |= {"I0", "ndiode"}
if "tfe" in all_modes:
used |= {"A_tfe", "E00"}
if "fn" in all_modes:
used |= {"A_fn", "B_fn"}
if "sclc" in all_modes:
used |= {"K_sclc"}
return used
[ドキュメント]
def junction_current_forward(Vd, p, mode_set, temperature):
"""
フォワードバイアス条件における総接合電流を計算します。
`mode_set` で指定されたメカニズムに基づいて、対応する電流成分（ダイオード、TFE、FN、SCLC）を合計します。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param p: パラメータ値を含む辞書。
:type p: dict
:param mode_set: フォワードバイアスで使用するメカニズムのセット。
:type mode_set: set
:param temperature: 温度 (K)。
:type temperature: float
:returns: 総接合電流。
:rtype: numpy.ndarray or float
"""
j = 0.0
if "diode" in mode_set:
j += j_diode(Vd, p["I0"], p["ndiode"], temperature)
if "tfe" in mode_set:
j += j_tfe_forward(Vd, p["A_tfe"], p["E00"])
if "fn" in mode_set:
j += j_fn_forward(Vd, p["A_fn"], p["B_fn"])
if "sclc" in mode_set:
j += j_sclc(Vd, p["K_sclc"])
return j
[ドキュメント]
def junction_current_reverse(Vd, p, mode_set, temperature):
"""
リバースバイアス条件における総接合電流を計算します。
`mode_set` で指定されたメカニズムに基づいて、対応する電流成分（ダイオード、TFE、FN、SCLC）を合計します。
:param Vd: ダイオード接合にかかる電圧。
:type Vd: numpy.ndarray or float
:param p: パラメータ値を含む辞書。
:type p: dict
:param mode_set: リバースバイアスで使用するメカニズムのセット。
:type mode_set: set
:param temperature: 温度 (K)。
:type temperature: float
:returns: 総接合電流。
:rtype: numpy.ndarray or float
"""
j = 0.0
if "diode" in mode_set:
j += j_diode(Vd, p["I0"], p["ndiode"], temperature)
if "tfe" in mode_set:
j += j_tfe_reverse(Vd, p["A_tfe"], p["E00"])
if "fn" in mode_set:
j += j_fn_reverse(Vd, p["A_fn"], p["B_fn"])
if "sclc" in mode_set:
j += j_sclc(Vd, p["K_sclc"])
return j
# =============================================================================
# model
# =============================================================================
[ドキュメント]
def same_sign(a, b):
"""
2つの数値が同じ符号を持つかどうかを判定します。
両方が正または両方が負の場合にTrueを返します。
片方または両方がゼロの場合の動作は、厳密には考慮されていません。
:param a: 1つ目の数値。
:type a: float
:param b: 2つ目の数値。
:type b: float
:returns: `a` と `b` が同じ符号を持つ場合はTrue、そうでない場合はFalse。
:rtype: bool
"""
return (a > 0.0 and b > 0.0) or (a < 0.0 and b < 0.0)
[ドキュメント]
def model(V, p, forwardIV="diode", reverseIV="diode+tfe", dV=0.05, temperature=300.0):
"""
拡張一ダイオードモデル（SDM）に基づいて電流-電圧（IV）カーブを計算します。
与えられた外部電圧 `V`、パラメータ `p`、および選択された輸送メカニズム (`forwardIV`, `reverseIV`) を使用して
対応する電流 `I` を計算します。接合電圧 `Vd` を未知数とする根探しアルゴリズム（Brentのメソッド）を
用いて安定的にモデル電流を求めます。
フォワードとリバースで異なるモデルが指定された場合、シグモイド関数 (`sigmoid_blend`) を使って
`Vd` が0V付近でスムーズにブレンドされます。
:param V: 外部印加電圧の配列。
:type V: numpy.ndarray
:param p: モデルパラメータを含む辞書。主要なキーは "I0", "ndiode", "IPV", "Rs", "Rsh", "A_tfe", "E00", "A_fn", "B_fn", "K_sclc"。
:type p: dict
:param forwardIV: フォワードバイアス領域で使用する輸送メカニズムを示す文字列。例: "diode", "diode+tfe+sclc"。
:type forwardIV: str
:param reverseIV: リバースバイアス領域で使用する輸送メカニズムを示す文字列。例: "diode+tfe", "diode+fn"。
:type reverseIV: str
:param dV: フォワードモデルとリバースモデルをブレンドする際の電圧幅。Vd=0V付近での遷移の鋭さを制御します。
:type dV: float
:param temperature: シミュレーションの温度 (K)。
:type temperature: float
:returns: 各外部電圧 `V` に対応する総電流の配列。
:rtype: numpy.ndarray
"""
fset = parse_mech_mode(forwardIV)
rset = parse_mech_mode(reverseIV)
same_model = (fset == rset)
V = np.atleast_1d(np.asarray(V, dtype=float))
Itotal = np.empty_like(V)
Rs = max(float(p["Rs"]), EPS_R)
Rsh = max(float(p["Rsh"]), EPS_R)
IPV = float(p["IPV"])
def residual_vd(vd_trial, v_target):
"""
接合電圧 `vd_trial` の残差関数を計算します。
この関数は `brentq` によって `vd_trial` の根を見つけるために使用されます。
"""
if same_model:
ij = junction_current_forward(vd_trial, p, fset, temperature)
else:
w = float(sigmoid_blend(vd_trial, dV))
ij_f = junction_current_forward(vd_trial, p, fset, temperature)
ij_r = junction_current_reverse(vd_trial, p, rset, temperature)
ij = w * ij_f + (1.0 - w) * ij_r
ibranch = ij + vd_trial / Rsh
v_calc = vd_trial + Rs * ibranch
return v_calc - v_target
for i, v_val in enumerate(V):
try:
# 根探しの初期探索範囲を設定。VdはVよりRs*I落ちるため、広めに設定。
width = max(5.0, abs(v_val) + Rs * (abs(IPV) + abs(v_val) / Rsh + 1.0))
lo = v_val - width
hi = v_val + width
# BrentQ のための探索範囲 [lo, hi] を確保。f(lo) と f(hi) が異なる符号を持つように調整。
flo = residual_vd(lo, v_val)
fhi = residual_vd(hi, v_val)
ntry = 0
while same_sign(flo, fhi) and ntry < 25:
width *= 2.0
lo = v_val - width
hi = v_val + width
flo = residual_vd(lo, v_val)
fhi = residual_vd(hi, v_val)
ntry += 1
if same_sign(flo, fhi): # 符号が異ならない場合、根が見つからなかったとみなし、近似値を返す
vd_root = v_val
else:
vd_root = brentq(
residual_vd, lo, hi,
args=(v_val,),
maxiter=300, xtol=1.0e-12, rtol=1.0e-10
)
# 根が見つかったVdから電流を再計算
if same_model:
ij = junction_current_forward(vd_root, p, fset, temperature)
else:
w = float(sigmoid_blend(vd_root, dV))
ij_f = junction_current_forward(vd_root, p, fset, temperature)
ij_r = junction_current_reverse(vd_root, p, rset, temperature)
ij = w * ij_f + (1.0 - w) * ij_r
ibranch = ij + vd_root / Rsh
Itotal[i] = ibranch - IPV
except Exception:
# 根探しに失敗した場合のフォールバック（例：V = 0 で Rs + Rsh がゼロに近い場合など）
Itotal[i] = v_val / max(Rs + Rsh, EPS_R) - IPV
return Itotal
# =============================================================================
# initial parameter estimation
# =============================================================================
[ドキュメント]
def estimate_voc_local_quadratic(V, I, nlsq_points=5):
"""
IVデータから開回路電圧（VOC）を局所二次多項式フィッティングで推定します。
電流がゼロに最も近い点の近傍データを使用して、二次多項式をフィッティングし、
その多項式の根からVOCを求めます。データ点が少ない場合は、線形補間を使用します。
:param V: 電圧データの配列。
:type V: numpy.ndarray
:param I: 電流データの配列。
:type I: numpy.ndarray
:param nlsq_points: 局所フィッティングに使用するデータ点の数。
:type nlsq_points: int
:returns: 推定された開回路電圧VOC。
:rtype: float
"""
V = np.asarray(V, dtype=float)
I = np.asarray(I, dtype=float)
if len(V) < 3:
return float(V[np.argmin(np.abs(I))])
nlsq_points = max(3, int(nlsq_points))
nlsq_points = min(nlsq_points, len(V))
ic = int(np.argmin(np.abs(I))) # abs(I) 最小点
half = nlsq_points // 2
i0 = max(0, ic - half)
i1 = min(len(V), i0 + nlsq_points)
i0 = max(0, i1 - nlsq_points)
x = V[i0:i1]
y = I[i0:i1]
deg = min(2, len(x) - 1)
coef = np.polyfit(x, y, deg)
if deg == 1:
a, b = coef
if abs(a) < 1.0e-30:
return float(V[ic])
return float(-b / a)
a, b, c = coef
if abs(a) < 1.0e-30:
if abs(b) < 1.0e-30:
return float(V[ic])
return float(-c / b)
roots = np.roots([a, b, c])
real_roots = roots[np.isreal(roots)].real
if len(real_roots) == 0:
return float(V[ic])
# 近傍窓の中心に最も近い根を採用
x_center = V[ic]
return float(real_roots[np.argmin(np.abs(real_roots - x_center))])
[ドキュメント]
def estimate_initial_params(V, I_meas, temperature):
"""
IVデータからSDM（拡張一ダイオードモデル）の初期パラメータを推定します。
データ平滑化、dI/dV解析、多項式根探し、負電圧領域でのTFE近似を用いて、
最適化に適した頑健な初期値を導出します。
:param V: 測定された電圧データ点の配列。
:type V: numpy.ndarray
:param I_meas: 測定された電流データ点の配列。
:type I_meas: numpy.ndarray
:param temperature: 測定時の温度 (K)。
:type temperature: float
:returns: 推定された初期パラメータを含む辞書。
:rtype: dict
"""
print("[Process] Estimating initial diode/TFE parameters from data...")
# 重複するV値を処理
uV, idx = np.unique(V, return_index=True)
uI = I_meas[idx]
# データ平滑化と微分
y_sm = smooth_polyfit(uI, window_points=5, poly_order=3)
dydv = np.gradient(y_sm, uV)
# Rsの推定: V>0 で dI/dV が最大となる点から
pos_mask = uV > 0.0
if np.any(pos_mask):
dpos = dydv[pos_mask]
xpos = uV[pos_mask]
idx_pos_local = int(np.argmax(dpos))
V_rs = float(xpos[idx_pos_local])
slope_pos = float(dpos[idx_pos_local])
else: # V>0 のデータがない場合、全体の最大勾配
V_rs = float(uV[-1])
slope_pos = float(np.max(dydv))
Rs_est = 1.0 / abs(slope_pos)
# Rshの推定: V<0 で dI/dV が最小（最も負）となる点から
neg_mask = uV < 0.0
if np.any(neg_mask):
dneg = dydv[neg_mask]
xneg = uV[neg_mask]
idx_neg_local = int(np.argmin(dneg))
V_sh = float(xneg[idx_neg_local])
slope_neg = float(dneg[idx_neg_local])
else: # V<0 のデータがない場合、全体の最小勾配
V_sh = float(uV[0])
slope_neg = float(np.min(dydv))
Rsh_est = 1.0 / max(abs(slope_neg), 1.0e-12)
# ISC, IPV, VOC の推定
ISC_est = local_poly_value(uV, y_sm, 0.0, npts=7, order=3)
IPV_est = max(-ISC_est, EPS_I)
VOC_est = solve_root_poly(uV, y_sm, target_y=0.0, order=3)
VOC_est = estimate_voc_local_quadratic(uV, y_sm, nlsq_points=5)
# Ish, I0 の推定
Ish_est = local_poly_value(uV, y_sm, V_sh, npts=7, order=3)
I0_est = max(ISC_est - Ish_est, 1.0e-12) # I0は通常正の値
# TFE 初期値推定 (V<0 の十分負側で ln(|I|)-V を一次近似)
Vtfe_est = np.nan
Itfe_est = np.nan
A_tfe_est = 1.0e-30
E00_est = 1.0e-3
neg_idx = np.where(uV < 0.0)[0]
if len(neg_idx) >= 5: # 負電圧領域に十分なデータ点がある場合
nneg = len(neg_idx)
ntail = max(5, int(np.ceil(0.4 * nneg))) # 最も負側の40%または5点以上を使用
tail_idx = neg_idx[:ntail] # 最も負側の領域
V_tail = uV[tail_idx]
I_tail = np.abs(y_sm[tail_idx])
mask_good = np.isfinite(V_tail) & np.isfinite(I_tail) & (I_tail > EPS_I)
V_tail = V_tail[mask_good]
I_tail = I_tail[mask_good]
if len(V_tail) >= 2: # 少なくとも2点あれば直線フィッティング可能
lnI = np.log(I_tail)
slope_tfe, intercept_tfe = np.polyfit(V_tail, lnI, 1)
if abs(slope_tfe) > 1.0e-12:
E00_est = 1.0 / abs(slope_tfe)
else:
E00_est = 1.0e-3 # デフォルト値
irep = int(np.argmin(V_tail)) # 最も負の電圧点
Vtfe_est = float(V_tail[irep])
Itfe_est = float(I_tail[irep])
# ln|I| = ln(A_tfe) + (-V)/E00 for reverse-TFE-like growth
# |I| = A_tfe * exp(-V/E00) より
A_tfe_est = Itfe_est * np.exp(Vtfe_est / max(E00_est, 1.0e-12))
A_tfe_est = max(float(A_tfe_est), 1.0e-30)
E00_est = max(float(E00_est), 1.0e-6)
# 推定されたパラメータとデフォルト値を結合
return {
"I0": I0_est,
"ndiode": 1.5,
"IPV": IPV_est,
"Rs": max(Rs_est, 1.0e-6), # Rsはゼロにならないようにする
"Rsh": max(Rsh_est, 1.0e3), # Rshは十分に大きい値から始める
"A_tfe": A_tfe_est,
"E00": E00_est,
"A_fn": 1.0e-30, # FN/SCLCは初期値推定でカバーされないため極小値
"B_fn": 10.0,
"K_sclc": 1.0e-30,
"VOC_est": VOC_est,
"ISC_est": ISC_est,
"Vsh_est": V_sh,
"Ish_est": Ish_est,
"Vtfe_est": Vtfe_est,
"Itfe_est": Itfe_est,
"dIdV_pos_rep": slope_pos,
"dIdV_neg_rep": slope_neg,
}
# =============================================================================
# fit support
# =============================================================================
[ドキュメント]
def pack_free_params(params, fix_set):
"""
辞書形式の全パラメータから、最適化対象の自由パラメータのみをNumPy配列にパックします。
`fix_set` に含まれるパラメータは除外されます。
`LOG_PARAMS` に含まれるパラメータは `log10` 変換されて配列に格納されます。
:param params: 全パラメータ（固定および自由）を含む辞書。
:type params: dict
:param fix_set: 固定されるパラメータ名のセット。
:type fix_set: set
:returns: 最適化に使用される自由パラメータのNumPy配列。
:rtype: numpy.ndarray
"""
p0 = []
for name in PARAM_NAMES:
if name in fix_set:
continue
val = float(params[name])
if name in LOG_PARAMS:
p0.append(np.log10(max(val, 1.0e-30))) # ゼロや負の値を避ける
else:
p0.append(val)
return np.asarray(p0, dtype=float)
[ドキュメント]
def unpack_free_params(p_free, base_params, fix_set):
"""
自由パラメータのNumPy配列を、固定パラメータと組み合わせて辞書形式の全パラメータにアンパックします。
`base_params` は固定パラメータの値を決定するために使用されます。
`LOG_PARAMS` に含まれるパラメータは `10**` 変換されて辞書に格納されます。
:param p_free: 最適化中の自由パラメータのNumPy配列。
:type p_free: numpy.ndarray
:param base_params: 固定パラメータの基準値を含む辞書。
:type base_params: dict
:param fix_set: 固定されたパラメータ名のセット。
:type fix_set: set
:returns: 更新された全パラメータを含む辞書。
:rtype: dict
"""
curr = dict(base_params) # base_params をコピーして、固定されていないパラメータを更新
idx = 0
for name in PARAM_NAMES:
if name in fix_set:
continue
val = float(p_free[idx])
if name in LOG_PARAMS:
curr[name] = 10.0 ** val
else:
curr[name] = val
idx += 1
return curr
[ドキュメント]
def signed_log10_current(I, eps=EPS_I):
"""
符号を保持したまま電流を圧縮する補助関数。
`sign(I) * log10(abs(I)+eps)` を返します。
電流値がゼロに近い場合に計算が不安定になることを防ぐため、`abs(I)` に `eps` を加えます。
:param I: 電流値の配列または単一値。
:type I: numpy.ndarray or float
:param eps: 絶対値計算時に加算される微小量。ゼロ除算や対数関数の引数が負になるのを防ぐ。
:type eps: float
:returns: 符号を保持したログスケールの電流値。
:rtype: numpy.ndarray or float
"""
I = np.asarray(I, dtype=float)
return np.sign(I) * np.log10(np.abs(I) + eps)
[ドキュメント]
def objective(p_free, V, I_meas, args, base_params, fix_set):
"""
最適化のための目的関数（残差平方和）を計算します。
測定電流とモデル電流の符号を保持した圧縮表現
`sign(I) * log10(abs(I)+EPS_I)` の残差平方和を使用します。
:param p_free: 最適化中の自由パラメータのNumPy配列。
:type p_free: numpy.ndarray
:param V: 測定された電圧データ点の配列。
:type V: numpy.ndarray
:param I_meas: 測定された電流データ点の配列。
:type I_meas: numpy.ndarray
:param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト（例: `forwardIV`, `reverseIV`）。
:type args: argparse.Namespace
:param base_params: 固定パラメータの基準値を含む辞書。
:type base_params: dict
:param fix_set: 固定されたパラメータ名のセット。
:type fix_set: set
:returns: 目的関数の値（残差平方和）。
:rtype: float
"""
curr = unpack_free_params(p_free, base_params, fix_set)
I_calc = model(
V, curr,
forwardIV=args.forwardIV,
reverseIV=args.reverseIV,
dV=args.dv_blend,
temperature=args.temperature,
)
# 符号付き log 圧縮スケールでの残差平方和を計算。
y_meas = signed_log10_current(I_meas, EPS_I)
y_calc = signed_log10_current(I_calc, EPS_I)
return float(np.sum((y_meas - y_calc) ** 2))
# =============================================================================
# error estimation
# =============================================================================
[ドキュメント]
def get_jacobian(p_free, V, temperature, args, base_params, fix_set):
"""
モデルのヤコビ行列を数値的に計算します。
各自由パラメータを微小量 `eps` だけ摂動させ、`log10(|I|)` の変化を計算して
偏微分（ヤコビ行列の要素）を近似します。
:param p_free: 最適化された自由パラメータのNumPy配列。
:type p_free: numpy.ndarray
:param V: 測定された電圧データ点の配列。
:type V: numpy.ndarray
:param temperature: シミュレーションの温度 (K)。
:type temperature: float
:param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト。
:type args: argparse.Namespace
:param base_params: 固定パラメータの基準値を含む辞書。
:type base_params: dict
:param fix_set: 固定されたパラメータ名のセット。
:type fix_set: set
:returns: 各データ点に対する各自由パラメータの偏微分からなるヤコビ行列。
:rtype: numpy.ndarray
"""
eps = 1.0e-5 # 微小摂動量
n_params = len(p_free)
n_points = len(V)
J = np.zeros((n_points, n_params), dtype=float)
def get_log_model(pv):
"""モデル電流の符号付きlog圧縮表現を計算するヘルパー関数。"""
curr = unpack_free_params(pv, base_params, fix_set)
I_calc = model(
V, curr,
forwardIV=args.forwardIV,
reverseIV=args.reverseIV,
dV=args.dv_blend,
temperature=temperature,
)
return signed_log10_current(I_calc, EPS_I)
f0 = get_log_model(p_free) # 現在のパラメータでのモデル出力
# 各自由パラメータに対して数値微分を実行
for i in range(n_params):
pp = np.copy(p_free)
pp[i] += eps # i番目のパラメータを微小量だけ変更
fi = get_log_model(pp)
J[:, i] = (fi - f0) / eps # 偏微分を近似
return J
[ドキュメント]
def estimate_errors(res, V, I_meas, temperature, args, base_params, fix_set):
"""
フィッティング結果からパラメータの誤差とモデル電流の信頼区間を推定します。
ヤコビ行列と残差平方和を用いて共分散行列を計算し、誤差伝播によって
各パラメータの誤差とモデル電流の `log10` スケールでの標準偏差を算出します。
:param res: `scipy.optimize.minimize` の結果オブジェクト。最適化されたパラメータ (`res.x`) と目的関数の値 (`res.fun`) を含む。
:type res: scipy.optimize.OptimizeResult
:param V: 測定された電圧データ点の配列。
:type V: numpy.ndarray
:param I_meas: 測定された電流データ点の配列。
:type I_meas: numpy.ndarray
:param temperature: シミュレーションの温度 (K)。
:type temperature: float
:param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト。
:type args: argparse.Namespace
:param base_params: 固定パラメータの基準値を含む辞書。
:type base_params: dict
:param fix_set: 固定されたパラメータ名のセット。
:type fix_set: set
:returns: errors_dict (各パラメータ名とその推定誤差の辞書), sigma_log_I (モデル電流のlog10値の標準偏差の配列)
:rtype: tuple[dict, numpy.ndarray]
"""
p_opt = np.asarray(res.x, dtype=float) # 最適化された自由パラメータ
n = len(I_meas) # データ点数
p = len(p_opt) # 自由パラメータ数
if p == 0 or n <= p: # 自由パラメータがない、またはデータ点数が少なすぎる場合
return {}, np.zeros(len(V), dtype=float)
J = get_jacobian(p_opt, V, temperature, args, base_params, fix_set) # ヤコビ行列
sigma2_hat = float(res.fun) / max(n - p, 1) # 残差平方和から残差分散を推定
try:
# 共分散行列の計算: (J^T J)^-1 * sigma_hat^2
# np.linalg.pinv は特異行列の場合も対応できる擬似逆行列
cov = sigma2_hat * np.linalg.pinv(J.T @ J)
# 自由パラメータ自体の誤差（logスケールの場合）
p_errors_internal = np.sqrt(np.maximum(np.diag(cov), 0.0))
# モデル電流のlog10値の標準偏差
sigma_log_I = np.sqrt(np.maximum(np.diag(J @ cov @ J.T), 0.0))
errors_dict = {}
idx = 0
for name in PARAM_NAMES:
if name in fix_set:
continue
if name in LOG_PARAMS:
# log10変換されたパラメータの誤差を元のスケールに戻す
val = 10.0 ** p_opt[idx]
errors_dict[name] = val * np.log(10.0) * p_errors_internal[idx]
else:
errors_dict[name] = p_errors_internal[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), dtype=float)
# =============================================================================
# plotting
# =============================================================================
[ドキュメント]
def initialize_plot(V, I_meas, I_sim, title = None):
"""
最適化プロセスのための初期プロットを設定し、アニメーションオブジェクトを返します。
リアルタイムでのフィッティング進捗状況を表示するために使用されます。
線形スケールと対数絶対値スケールの2つのサブプロットが作成されます。
:param V: 電圧データの配列。
:type V: numpy.ndarray
:param I_meas: 測定された電流データの配列。
:type I_meas: numpy.ndarray
:param I_sim: 初期シミュレーション電流データの配列。
:type I_sim: numpy.ndarray
:param title: プロットのメインタイトル。Noneの場合、タイトルは設定されない。
:type title: str or None
:returns: fig (matplotlib.figure.Figure), axes (matplotlib.axes.Axesの配列), anim (アニメーションに使うオブジェクトの辞書)
:rtype: tuple[matplotlib.figure.Figure, numpy.ndarray, dict]
"""
anim = {}
plt.ion() # インタラクティブモードをオン
fig, axes = plt.subplots(1, 2, figsize=(11, 5))
ax1, ax2 = axes
if title: ax1.set_title(title)
ax1.plot(V, I_meas, "o", color="black", alpha=0.3, markersize=4, label="data")
line0, = ax1.plot(V, I_sim, "-", color="blue", lw=0.5, label="initial")
line1, = ax1.plot(V, I_sim, "-", color="red", lw=1, label="sim")
band1 = ax1.fill_between(V, I_sim, I_sim, color="lightblue", alpha=0.5) # 初期状態では幅ゼロ
ax1.axhline(0.0, color="k", lw=0.5)
ax1.axvline(0.0, color="k", lw=0.5)
ax1.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
ax1.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
ax1.set_xlabel("V / V", fontsize = fontsize)
ax1.set_ylabel("I / A", fontsize = fontsize)
ax1.grid(True)
ax1.legend()
ax2.set_title("Log |I|")
ax2.semilogy(V, np.abs(I_meas) + EPS_I, "o", color="black", alpha=0.3, markersize=4)
line02, = ax2.semilogy(V, np.abs(I_sim) + EPS_I, "-", color="blue", lw=0.5)
line2, = ax2.semilogy(V, np.abs(I_sim) + EPS_I, "-", color="red", lw=1)
band2 = ax2.fill_between(V, np.abs(I_sim) + EPS_I, np.abs(I_sim) + EPS_I, color="lightblue", alpha=0.5) # 初期状態では幅ゼロ
ax2.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
ax2.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
ax2.set_xlabel("V / V", fontsize = fontsize)
ax2.set_ylabel("I / A", fontsize = fontsize)
ax2.grid(True, which="both")
plt.tight_layout()
anim = {
"fig": fig,
"ax1": ax1,
"ax2": ax2,
"line1": line1,
"line2": line2,
"band1": band1,
"band2": band2,
}
return fig, axes, anim
[ドキュメント]
def plot_iv(fig, axes, V, I_meas, I_start, I_final, info, mode, sigma_log=None, mode_label = "", outfile=None, pause=False):
"""
測定データとモデルIVカーブをプロットします。
線形スケールと対数絶対値スケールの2つのサブプロットを作成します。
`sigma_log` が提供された場合、モデル電流の信頼区間を水色領域で表示します。
:param fig: プロットに使うmatplotlib.figure.Figureオブジェクト。Noneの場合、新しいFigureを作成。
:type fig: matplotlib.figure.Figure or None
:param axes: プロットに使うmatplotlib.axes.Axesオブジェクトの配列。Noneの場合、新しいAxesを作成。
:type axes: numpy.ndarray or None
:param V: プロットする電圧データの配列。
:type V: numpy.ndarray
:param I_meas: 測定された電流データの配列。Noneの場合、測定データはプロットされない。
:type I_meas: numpy.ndarray or None
:param I_start: 初期モデルによって計算された電流データの配列。Noneの場合、初期モデルはプロットされない。
:type I_start: numpy.ndarray or None
:param I_final: 最終モデルによって計算された電流データの配列。
:type I_final: numpy.ndarray
:param info: ファイル名やレコード時間などの情報を含む辞書。プロットのタイトルに使用。
:type info: dict
:param mode: 現在の実行モード（例: "init", "fit", "sim"）を示す文字列。プロットのタイトルに使用。
:type mode: str
:param sigma_log: モデル電流の `log10` 値の標準偏差の配列。信頼区間を表示するために使用。Noneの場合、信頼区間は表示されない。
:type sigma_log: numpy.ndarray or None
:param mode_label: プロットのタイトルに追加するモードのラベル。
:type mode_label: str
:param outfile: プロットを保存するファイルパス。Noneの場合、保存しない。
:type outfile: str or None
:returns: None
:rtype: None
"""
if fig is None:
fig, axes = plt.subplots(1, 2, figsize=(11, 5))
ax = axes[0]
ax.set_title(f"{info['FileName']} ({mode_label})")
if I_meas is not None:
ax.plot(V, I_meas, "o", color="black", alpha=0.3, markersize=4, label="data")
if I_start is not None:
ax.plot(V, I_start, "-", color="blue", lw=0.5, label="initial")
ax.plot(V, I_final, "-", color="red", lw=1, label="final")
if sigma_log is not None:
log_abs = np.log10(np.abs(I_final) + EPS_I)
upper = 10.0 ** (log_abs + sigma_log)
lower = 10.0 ** (log_abs - sigma_log)
sign = np.sign(I_final) # 電流の符号を保持
ax.fill_between(V, sign * lower, sign * upper, color="lightblue", alpha=0.5)
ax.axhline(0.0, color="k", lw=0.5)
ax.axvline(0.0, color="k", lw=0.5)
ax.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
ax.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
ax.set_xlabel("V / V", fontsize = fontsize)
ax.set_ylabel("I / A", fontsize = fontsize)
ax.grid(True)
ax.legend()
ax = axes[1]
ax.set_title("Log |I|")
if I_meas is not None:
ax.semilogy(V, np.abs(I_meas) + EPS_I, "o", color="black", alpha=0.3, markersize=4)
ax.semilogy(V, np.abs(I_final) + EPS_I, "-", color="red", lw=1)
if sigma_log is not None:
log_abs = np.log10(np.abs(I_final) + EPS_I)
upper = 10.0 ** (log_abs + sigma_log)
lower = 10.0 ** (log_abs - sigma_log)
ax.fill_between(V, lower, upper, color="lightblue", alpha=0.5)
ax.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
ax.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
ax.set_xlabel("V / V", fontsize = fontsize)
ax.set_ylabel("I / A", fontsize = fontsize)
ax.grid(True, which="both")
fig.tight_layout()
if outfile:
print(f"[I/O] Save plot: {outfile}")
fig.savefig(outfile, dpi=160)
if pause:
plt.show()
else:
plt.pause(0.01)
[ドキュメント]
def setup_parameters(params, args, csv_fix, args_fix, forwardIV, reverseIV):
"""
パラメータ辞書を最終的に設定し、固定されるパラメータを決定します。
コマンドライン引数で上書きされた値、CSVファイルから読み込まれた値、
デフォルト値、および未使用のモデルメカニズムに基づいて自動的に固定されるパラメータを統合します。
:param params: 現在のパラメータ値を含む辞書。
:type params: dict
:param args: 解析されたコマンドライン引数を含むオブジェクト。
:type args: argparse.Namespace
:param csv_fix: CSVファイルで固定されたパラメータ名のセット。
:type csv_fix: set
:param args_fix: コマンドライン引数で固定されたパラメータ名のリスト。
:type args_fix: list
:param forwardIV: フォワードバイアスモデルのメカニズムを示す文字列。
:type forwardIV: str
:param reverseIV: リバースバイアスモデルのメカニズムを示す文字列。
:type reverseIV: str
:returns: params (最終的なパラメータ値の辞書), fix_set (最終的に固定されるパラメータ名のセット), used (使用されるパラメータ名のセット), auto_fixed (自動的に固定されたパラメータ名のセット)
:rtype: tuple[dict, set, set, set]
:raises ValueError: 不明な固定パラメータ名が指定された場合。
"""
# コマンドライン引数によるパラメータ上書き
for name in PARAM_NAMES:
val = getattr(args, name)
if val is not None:
params[name] = val
# 不足しているパラメータにデフォルト値を設定
defaults = {
"I0": 1.0e-10,
"ndiode": 1.5,
"IPV": 1.0e-6,
"Rs": 1.0,
"Rsh": 1.0e9,
"A_tfe": 1.0e-30,
"E00": 1.0e-3,
"A_fn": 1.0e-30,
"B_fn": 10.0,
"K_sclc": 1.0e-30,
}
for k, v in defaults.items():
params.setdefault(k, v)
# 固定パラメータの決定
fix_set = set(csv_fix) | set(args_fix) # CSVとコマンドライン引数の両方で指定されたものを結合
# 使用されていないモデルパラメータを自動的に固定
used = used_params_from_modes(forwardIV, reverseIV)
auto_fixed = set(PARAM_NAMES) - used
fix_set |= auto_fixed
# 不明な固定パラメータ名のチェック
unknown_fix = [x for x in fix_set if x not in PARAM_NAMES]
if unknown_fix:
raise ValueError(f"Unknown parameter name(s) in --fix: {unknown_fix}")
print(f"Used parameters : {' '.join(sorted(used))}")
print(f"Auto-fixed unused : {' '.join(sorted(auto_fixed)) if auto_fixed else '(none)'}")
return params, fix_set, used, auto_fixed
[ドキュメント]
def print_parameter_repeat(params, fix_set, used, auto_fixed, title=""):
"""
現在のモデルパラメータとその固定状態をコンソールに出力します。
これは、`print_args_and_derived` の「Derived values」セクションに似ていますが、
最終的に決定された各パラメータの具体的な値と、それが固定されているかどうかの状態を詳細に示します。
:param params: 最終的なパラメータ値を含む辞書。
:type params: dict
:param fix_set: 固定されるパラメータ名のセット。
:type fix_set: set
:param used: モデルで使用されるパラメータ名のセット。
:type used: set
:param auto_fixed: モデルによって自動的に固定されたパラメータ名のセット。
:type auto_fixed: set
:param title: 出力の前に表示されるタイトル文字列。
:type title: str
:returns: None
:rtype: None
"""
if title:
print(title)
for name in PARAM_NAMES:
status = ""
if name in fix_set:
status = "(Fixed)"
if name in auto_fixed:
status = "(Auto-Fixed)"
print(f" {name:8s}: {params.get(name, 'N/A'):.6e} {status}")
print()
# =============================================================================
# fit execution
# =============================================================================
[ドキュメント]
def exec_fit(V, I_meas, info, csv_path, params, fix_set, args):
"""
fit モードの本体処理。
最適化、誤差推定、保存、最終プロットまでをまとめて実行する。
この関数は、`scipy.optimize.minimize` を用いてSDMパラメータを測定データにフィッティングします。
フィッティング中には、コンソールへの進捗表示とリアルタイムのIVカーブアニメーションが行われます。
最適化後、パラメータの誤差が推定され、結果がCSVファイルに保存され、最終的なプロットが表示されます。
:param V: 測定された電圧データ点の配列。
:type V: numpy.ndarray
:param I_meas: 測定された電流データ点の配列。
:type I_meas: numpy.ndarray
:param info: ファイル名やレコード時間などの情報を含む辞書。プロットのタイトルに使用。
:type info: dict
:param csv_path: パラメータを保存するCSVファイルのパス。
:type csv_path: pathlib.Path
:param params: 初期パラメータ値を含む辞書。
:type params: dict
:param fix_set: 固定されるパラメータ名のセット。
:type fix_set: set
:param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト。
:type args: argparse.Namespace
:returns: params (最適化後のパラメータ値の辞書), rss (残差平方和), errors (各パラメータの推定誤差の辞書), sigma_log_curve (モデル電流のlog10値の標準偏差の配列)
:rtype: tuple[dict, float, dict, numpy.ndarray]
"""
print(f"Execute {args.mode} mode.")
params_init = dict(params)
I_start = model(V, params, args.forwardIV, args.reverseIV, args.dv_blend, args.temperature)
fig, axes, anim = initialize_plot(V, I_meas, I_start, title="Optimization progress")
p0 = pack_free_params(params, fix_set)
state = {"count": 0}
def callback(xk):
state["count"] += 1
curr = unpack_free_params(xk, params, fix_set)
if state["count"] % max(args.ninterval_print, 1) == 0:
rss = objective(xk, V, I_meas, args, params, fix_set)
pstr = " ".join(
f"{n}={curr[n]:.3e}" for n in PARAM_NAMES if n not in fix_set
)
print(f"{state['count']:4d}: RSS={rss:.6e} {pstr}")
if state["count"] % max(args.ninterval_plot, 1) == 0:
I_anim = model(
V, curr,
forwardIV=args.forwardIV,
reverseIV=args.reverseIV,
dV=args.dv_blend,
temperature=args.temperature,
)
anim["line1"].set_ydata(I_anim)
anim["line2"].set_ydata(np.abs(I_anim) + EPS_I)
anim["ax1"].relim()
anim["ax1"].autoscale_view()
anim["ax2"].relim()
anim["ax2"].autoscale_view()
anim["fig"].canvas.draw()
anim["fig"].canvas.flush_events()
plt.pause(0.001)
if len(p0) == 0:
print("All parameters are fixed. Skipping optimization.")
class DummyResult:
x = np.array([], dtype=float)
fun = objective(np.array([], dtype=float), V, I_meas, args, params, fix_set)
res = DummyResult()
rss = float(res.fun)
else:
res = minimize(
objective,
p0,
args=(V, I_meas, args, params, fix_set),
method=args.method,
callback=callback,
options={
"maxiter": int(args.nmaxiter),
"xatol": float(args.tol),
"fatol": float(args.tol),
"disp": False,
},
)
params = unpack_free_params(res.x, params, fix_set)
rss = float(res.fun)
errors, sigma_log_curve = estimate_errors(
res, V, I_meas, args.temperature, args, params, fix_set
)
save_param_csv(csv_path, params, fix_set, errors=errors, rss=rss)
plt.ioff()
print("\n--- Fit Results ---")
print(f"RSS_logy : {rss:.6e}")
for n in PARAM_NAMES:
if n in errors:
print(f" {n:8s}: {params[n]:.6e} +/- {errors[n]:.2e}")
else:
print(f" {n:8s}: {params[n]:.6e}" + (" (Fixed)" if n in fix_set else ""))
print()
print("Parameters estimated from input IV data (Use if the fitting accurary is poor)")
est = estimate_initial_params(V, I_meas, args.temperature)
print(f" Rs_est : {est['Rs']:.6f} ohm")
print(f" Rsh_est : {est['Rsh']:.6f} ohm")
print(f" ISC_est : {est['ISC_est']:.6e} A")
print(f" VOC_est : {est['VOC_est']:.6f} V")
I_final = model(
V, params,
forwardIV=args.forwardIV,
reverseIV=args.reverseIV,
dV=args.dv_blend,
temperature=args.temperature,
)
# Excel保存
outfile_xlsx = f"{args.outprefix}-fitted.xlsx"
save_iv_to_excel(
outfile_xlsx,
V=V,
I_input=I_meas,
I_init=I_start,
I_fit=I_final,
info=info,
params_init=params_init,
params_final=params,
rss=rss,
errors=errors,
)
# 画像保存してから表示
outfile_iv = f"{args.outprefix}-iv.png"
plot_iv(
fig, axes, V, I_meas, I_start, I_final, info, args.mode,
sigma_log=sigma_log_curve,
outfile=outfile_iv,
pause=False
)
input("\nPress ENTER to terminate\n")
return params, rss, errors, sigma_log_curve
# =============================================================================
# main
# =============================================================================
_original_print = None
_redirect_fp = None
[ドキュメント]
def dual_print(*args, **kwargs):
_original_print(*args, **kwargs)
file_kwargs = dict(kwargs)
file_kwargs.pop("file", None) # file指定があっても無視してログファイルへ
file_kwargs.setdefault("flush", True)
_original_print(*args, file=_redirect_fp, **file_kwargs)
[ドキュメント]
def main():
"""
スクリプトのエントリーポイント。コマンドライン引数を解析し、IVフィッティングツールを実行します。
この関数は、初期値推定 (`init`)、モデルフィッティング (`fit`)、またはシミュレーション (`sim`) の
いずれかのモードで動作します。データ読み込み、パラメータ管理、最適化、結果のプロットを行います。
:returns: None
:rtype: None
"""
global _original_print, _redirect_fp
print()
print("Read commandline arguments")
args, parser = initialize()
validate_args(args)
_original_print = builtins.print
logfile = f"{args.outprefix}_{args.mode}.txt"
print(f"Open log file [{logfile}]")
_redirect_fp = open(logfile, "w", encoding="utf-8")
builtins.print = dual_print
print_args_and_derived(args)
print(f"forwardIV model : {args.forwardIV}")
print(f"reverseIV model : {args.reverseIV}")
# データ入力処理
print()
if args.infile:
print(f"Read data from [{args.infile}]")
print(f" x range: {args.xmin} - {args.xmax}, every {args.ndataskip} data")
V, I_meas, info = read_data(args.infile, xmin=args.xmin, xmax=args.xmax, ndataskip=args.ndataskip)
csv_path = Path(args.infile).stem + "-parameters.csv"
else:
if args.mode != "sim":
raise ValueError("--infile is required for mode init/fit/analyze.")
print(f"Create simulation data")
# シミュレーションモードの場合、ダミーデータを生成
V = np.linspace(-0.2, 0.8, 500)
I_meas = None
info = {"FileName": "Simulation", "RecordTime": "N/A"}
csv_path = "sim-parameters.csv"
# パラメータの読み込みと初期値設定
print()
if os.path.exists(csv_path):
print(f"Read parameters from [{csv_path}]")
params, csv_fix = load_param_csv(csv_path)
else:
print(f"Parameter file [{csv_path}] does not exist. Use default.")
params = {}
csv_fix = set()
if (args.mode == "init" or args.mode == "analyze") or not any(name in params for name in PARAM_NAMES):
if I_meas is not None:
print()
print(f"Estimate initial parameters at T={args.temperature}")
est = estimate_initial_params(V, I_meas, args.temperature)
else:
print()
print(f"Use default initial parameters at T={args.temperature}")
est = DEFAULT_PARAMETERS
params.update(est)
print()
print("Final setup model parameters from paramter file and commandline arguments")
params, fix_set, used, auto_fixed = setup_parameters(params, args, csv_fix, args.fix, args.forwardIV, args.reverseIV)
print_parameter_repeat(params, fix_set, used, auto_fixed, title="=== Final model parameters ===")
print()
if args.mode == "analyze" or args.mode == "init":
print(f"Execute {args.mode} mode.")
print(" Extracted parameters")
print(f" ISC_est : {params.get('ISC_est', np.nan):.6e} A")
print(f" VOC_est : {params.get('VOC_est', np.nan):.6f} V")
print(f" Vsh_est : {params.get('Vsh_est', np.nan):.6f} V")
print(f" Ish_est : {params.get('Ish_est', np.nan):.6e} A")
print(f" Vtfe_est: {params.get('Vtfe_est', np.nan):.6f} V")
print(f" Itfe_est: {params.get('Itfe_est', np.nan):.6e} A")
I_model = model(
V, params,
forwardIV=args.forwardIV,
reverseIV=args.reverseIV,
dV=args.dv_blend,
temperature=args.temperature,
)
if args.mode == "init":
save_param_csv(csv_path, params, fix_set)
outfile_xlsx = f"{args.outprefix}-fitted.xlsx"
save_iv_to_excel(
outfile_xlsx,
V=V,
I_input=I_meas,
I_init=None,
I_fit=I_model,
info=info,
params_init=params,
params_final=None,
rss=None,
errors=None,
)
plot_iv(
None, None, V, I_meas, None, I_model, info, args.mode,
sigma_log=None,
outfile=f"{args.outprefix}-iv.png",
pause=False,
)
input("\nPress ENTER to terminate\n")
elif args.mode == "fit":
params, rss, errors, sigma_log_curve = exec_fit(
V, I_meas, info, csv_path, params, fix_set, args
)
elif args.mode == "sim":
print(f"Execute {args.mode} mode.")
I_model = model(
V, params,
forwardIV=args.forwardIV,
reverseIV=args.reverseIV,
dV=args.dv_blend,
temperature=args.temperature,
)
outfile_xlsx = f"{args.outprefix}-fitted.xlsx"
save_iv_to_excel(
outfile_xlsx,
V=V,
I_input=None,
I_init=None,
I_fit=I_model,
info=info,
params_init=params,
params_final=None,
rss=None,
errors=None,
)
plot_iv(
None, None, V, I_meas, None, I_model, info, args.mode,
sigma_log=None,
outfile=f"{args.outprefix}-iv.png",
pause=False,
)
input("\nPress ENTER to terminate\n")
if __name__ == "__main__":
try:
main()
except Exception:
traceback.print_exc()
sys.exit(1)
finally:
if _redirect_fp: _redirect_fp.close()