pvfit.py ダウンロード/コピー
pvfit.py
pvfit.py
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4概要: 拡張一ダイオードモデル（SDM）IVフィッティングツール
5
6詳細説明:
7 このモジュールは、太陽電池や他の半導体デバイスの電流-電圧（IV）特性を、
8 拡張一ダイオードモデル（SDM）を用いてフィッティング、初期値推定、またはシミュレーションするためのツールを提供します。
9 複数の輸送メカニズム（ダイオード、トンネル誘起フォワードエミッション (TFE)、
10 非理想フォワード電流 (FN)、空間電荷制限電流 (SCLC)）を組み合わせて利用できます。
11
12主な機能:
13 フォワードバイアスとリバースバイアスのモデルを独立に選択可能。
14 ダイオード、TFE、FN、SCLCの各メカニズムを任意に組み合わせ可能。
15 数値解法は接合電圧 Vd を未知数とする根探しで安定化。
16 モデルに未使用のパラメータは自動的に固定。
17 forwardIV と reverseIV が同一の場合、再計算を回避して効率化。
18 最適化途中のIVカーブをアニメーション表示（常に有効）。
19 線形近似に基づくパラメータ誤差推定。
20 モデル電流の標準偏差から信頼性区間を水色領域で表示。
21
22mode=init の初期値推定の概要:
23 入力IVを電圧昇順にソート。
24 I-Vを局所多項式で平滑化。
25 dI/dV から Rs, Rsh を推定。
26 ISC = I(V=0) を V=0 近傍の多項式で推定。
27 IPV = -ISC。
28 VOC = I(V)=0 の多項式根で推定。
29 Ish = V<0 の代表点 Vsh における電流。
30 I0 = ISC - Ish （シャント電流は差し引かない）。
31 TFE 初期値: V<0 の十分負側で ln(abs(I))-V を一次近似し、E00 = 1/abs(slope)、
32 A_tfe = abs(I) * exp(+V/E00) を代表点から推定。
33 ndiode は固定初期値 1.5。
34 FN / SCLC の初期値は結果に影響しない極小値。
35 mode=init の初期値は fit 用の頑健な初期値です。
36 物理的に厳密な意味づけは mode=fit で精密化してください。
37
38関連リンク:
39 pvfit_usage
40"""
41
42import os
43import sys
44import argparse
45import builtins
46from pathlib import Path
47import traceback
48import csv
49from openpyxl import Workbook
50from openpyxl.styles import Font, PatternFill
51
52import matplotlib.pyplot as plt
53import numpy as np
54from scipy.optimize import brentq, minimize
55
56KB = 1.380649e-23
57E_CHARGE = 1.602176634e-19
58
59PARAM_NAMES = [
60 "I0", "ndiode", "IPV", "Rs", "Rsh",
61 "A_tfe", "E00",
62 "A_fn", "B_fn",
63 "K_sclc",
64]
65LOG_PARAMS = {
66 "I0", "IPV", "Rs", "Rsh",
67 "A_tfe", "E00",
68 "A_fn", "B_fn",
69 "K_sclc",
70}
71
72EPS_I = 1.0e-15
73EPS_R = 1.0e-12
74
75DEFAULT_PARAMETERS = {
76 "I0": 1.0e-10,
77 "ndiode": 1.5,
78 "IPV": 1.0e-6,
79 "Rs": 1.0,
80 "Rsh": 1.0e9,
81 "A_tfe": 1.0e-30,
82 "E00": 1.0e-3,
83 "A_fn": 1.0e-30,
84 "B_fn": 10.0,
85 "K_sclc": 1.0e-30,
86 "VOC_est": 0.5,
87 "ISC_est": -1.0e-6,
88 "Vsh_est": -0.1,
89 "Ish_est": -2.0e-6,
90 "Vtfe_est": -0.2,
91 "Itfe_est": 1.0e-12,
92 }
93
94fontsize = 16
95
96
97def initialize():
98 """
99 概要: コマンドライン引数パーサーを初期化し、引数を解析します。
100
101 詳細説明:
102 各種実行モード (analyze, init, fit, sim)、最適化手法、
103 入力ファイル、温度、データクリッピング範囲、スキップ間隔、
104 フォワード/リバースIVモデル、ブレンド電圧幅、プロット/プリント間隔、
105 最大イテレーション数、許容誤差、固定パラメータなどの引数を定義します。
106
107 引数:
108 なし
109
110 戻り値:
111 :returns: args (解析された引数を含むオブジェクト), parser (ArgumentParserオブジェクト)
112 :rtype: tuple[argparse.Namespace, argparse.ArgumentParser]
113 """
114 parser = argparse.ArgumentParser(description="Extended SDM IV fitting tool with diode/TFE/FN/SCLC.")
115 parser.add_argument("--mode", choices=["analyze", "init", "fit", "sim"], default="fit",
116 help="Execution mode (init: initial parameter estimation, fit: fit parameters to data, sim: simulate IV curve)")
117 parser.add_argument("--method", default="Nelder-Mead",
118 help="scipy.optimize.minimize method (e.g., 'Nelder-Mead', 'Powell', 'BFGS', 'L-BFGS-B')")
119 parser.add_argument("--infile", help="Input IV CSV file path")
120 parser.add_argument("--temperature", type=float, default=300.0,
121 help="Temperature in K for the model")
122 parser.add_argument("--xmin", type=float, help="Minimum voltage (V) for data clipping, fitting, and plotting")
123 parser.add_argument("--xmax", type=float, help="Maximum voltage (V) for data clipping, fitting, and plotting")
124 parser.add_argument("--ndataskip", type=int, default=0, help="# of data to skip reading")
125 parser.add_argument("--forwardIV", default="diode",
126 help="Transport mechanisms for forward-bias side, e.g., 'diode' or 'diode+tfe+sclc'")
127 parser.add_argument("--reverseIV", default="diode+tfe",
128 help="Transport mechanisms for reverse-bias side, e.g., 'diode+tfe' or 'diode+fn'")
129 parser.add_argument("--dv_blend", type=float, default=0.05,
130 help="Sigmoid transition width (V) between reverse and forward models for junction voltage Vd")
131 parser.add_argument("--ninterval_print", type=int, default=10,
132 help="Console print interval during fitting iterations")
133 parser.add_argument("--ninterval_plot", type=int, default=10,
134 help="Animation update interval during fitting iterations")
135 parser.add_argument("--nmaxiter", type=int, default=1000,
136 help="Maximum optimizer iterations")
137 parser.add_argument("--tol", type=float, default=1.0e-7,
138 help="Optimizer tolerance (xatol and fatol)")
139 parser.add_argument("--fix", nargs="*", default=[],
140 help="Space-separated list of parameter names to fix during fitting (e.g., 'I0 ndiode')")
141 parser.add_argument("--outprefix", default="pvfit", help="Output file prefix")
142
143 for p in PARAM_NAMES:
144 parser.add_argument(f"--{p}", type=float, help=f"Override initial value for {p}")
145
146 args = parser.parse_args()
147 args.animate_fit = True # 最適化中のアニメーション表示は常に有効
148
149 return args, parser
150
151# =============================================================================
152# argument repeat / validation
153# =============================================================================
154
155def print_args_and_derived(args):
156 """
157 概要: コマンドライン引数とその派生値をコンソールに出力します。
158
159 詳細説明:
160 ユーザーが指定した入力値と、プログラム内部で決定される追加の設定（例：使用されるパラメータ、自動固定されるパラメータ）
161 を一覧表示し、設定の確認を容易にします。
162
163 引数:
164 :param args: 解析されたコマンドライン引数を含むオブジェクト。
165 :type args: argparse.Namespace
166
167 戻り値:
168 :returns: なし
169 :rtype: None
170 """
171 print("=== Arguments ===")
172 print(f"mode = {args.mode}")
173 print(f"infile = {args.infile}")
174 print(f"temperature = {args.temperature} K")
175 print(f"xmin = {args.xmin} V" if args.xmin is not None else "xmin = None")
176 print(f"xmax = {args.xmax} V" if args.xmax is not None else "xmax = None")
177 print(f"ndataskip = {args.ndataskip}")
178 print(f"forwardIV = {args.forwardIV}")
179 print(f"reverseIV = {args.reverseIV}")
180 print(f"dv_blend = {args.dv_blend} V")
181 print(f"ninterval_print = {args.ninterval_print}")
182 print(f"ninterval_plot = {args.ninterval_plot}")
183 print(f"nmaxiter = {args.nmaxiter}")
184 print(f"tol = {args.tol}")
185 print(f"fix = {' '.join(args.fix) if args.fix else '(none)'}")
186 print(f"outprefix = {args.outprefix}")
187
188 for name in PARAM_NAMES:
189 val = getattr(args, name, None)
190 print(f"{name:12s} = {val}" if val is not None else f"{name:12s} = None")
191 print()
192
193 print("=== Derived values ===")
194 used = used_params_from_modes(args.forwardIV, args.reverseIV)
195 auto_fixed = set(PARAM_NAMES) - used
196 print(f"used_params = {' '.join(sorted(used)) if used else '(none)'}")
197 print(f"auto_fixed = {' '.join(sorted(auto_fixed)) if auto_fixed else '(none)'}")
198 print()
199
200
201def validate_args(args):
202 """
203 概要: コマンドライン引数の有効性を検証します。
204
205 詳細説明:
206 モードごとの必須引数、数値引数の範囲、パラメータ名の整合性などをチェックし、
207 無効な引数が存在する場合は ValueError を発生させます。
208
209 引数:
210 :param args: 解析されたコマンドライン引数を含むオブジェクト。
211 :type args: argparse.Namespace
212
213 戻り値:
214 :returns: なし
215 :rtype: None
216
217 例外:
218 :raises ValueError: 無効な引数が見つかった場合。
219 """
220 if args.mode in ("init", "fit", "analyze"):
221 if not args.infile:
222 raise ValueError(f"--infile is required for mode={args.mode}")
223
224 if args.ndataskip < 0:
225 raise ValueError(f"ndataskip must be >= 0, got {args.ndataskip}")
226
227 if args.xmin is not None and args.xmax is not None and args.xmin > args.xmax:
228 raise ValueError(f"xmin must be <= xmax, got xmin={args.xmin}, xmax={args.xmax}")
229
230 if args.temperature <= 0.0:
231 raise ValueError(f"temperature must be > 0 K, got {args.temperature}")
232
233 if args.dv_blend < 0.0:
234 raise ValueError(f"dv_blend must be >= 0 V, got {args.dv_blend}")
235
236 if args.ninterval_print <= 0:
237 raise ValueError(f"ninterval_print must be > 0, got {args.ninterval_print}")
238
239 if args.ninterval_plot <= 0:
240 raise ValueError(f"ninterval_plot must be > 0, got {args.ninterval_plot}")
241
242 if args.nmaxiter <= 0:
243 raise ValueError(f"nmaxiter must be > 0, got {args.nmaxiter}")
244
245 if args.tol <= 0.0:
246 raise ValueError(f"tol must be > 0, got {args.tol}")
247
248 unknown_fix = [x for x in args.fix if x not in PARAM_NAMES]
249 if unknown_fix:
250 raise ValueError(f"Unknown parameter name(s) in --fix: {unknown_fix}")
251
252 for name in LOG_PARAMS:
253 val = getattr(args, name, None)
254 if val is not None and val <= 0.0:
255 raise ValueError(f"{name} must be > 0 because it is log-scaled, got {val}")
256
257 # 非logパラメータ側
258 if getattr(args, "ndiode", None) is not None and args.ndiode <= 0.0:
259 raise ValueError(f"ndiode must be > 0, got {args.ndiode}")
260
261 if getattr(args, "B_fn", None) is not None and args.B_fn <= 0.0:
262 raise ValueError(f"B_fn must be > 0, got {args.B_fn}")
263
264
265# =============================================================================
266# data I/O
267# =============================================================================
268
269def read_data(infile, xmin=None, xmax=None, ndataskip=0):
270 """
271 概要: 指定されたCSVファイルからIVデータ（電圧と電流）を読み込みます。
272
273 詳細説明:
274 ファイルからメタデータを抽出し、電圧昇順にデータをソートし、
275 指定された電圧範囲でクリッピングし、データスキップ間隔を適用します。
276
277 引数:
278 :param infile: 入力IV CSVファイルのパス。
279 :type infile: str
280 :param xmin: データをクリップする最小電圧 (V)。Noneの場合、下限なし。
281 :type xmin: float or None
282 :param xmax: データをクリップする最大電圧 (V)。Noneの場合、上限なし。
283 :type xmax: float or None
284 :param ndataskip: 読み込み時にスキップするデータ点の数（n+1点ごとにデータを保持）。
285 :type ndataskip: int
286
287 戻り値:
288 :returns: V (読み込まれた電圧データのNumPy配列), I (読み込まれた電流データのNumPy配列), info (ファイル名と記録時間を含む辞書)
289 :rtype: tuple[numpy.ndarray, numpy.ndarray, dict]
290
291 例外:
292 :raises ValueError: 有効なデータが見つからない場合、またはクリッピング後にデータ点が残らない場合。
293 """
294 data_started = False
295 raw_v = []
296 raw_i = []
297 info = {"FileName": infile, "RecordTime": "Unknown"}
298
299 with open(infile, "r", encoding="utf-8-sig", newline="") as f:
300 reader = csv.reader(f)
301 for row in reader:
302 if not row:
303 continue
304
305 cols = [c.strip() for c in row]
306
307 if len(cols) >= 3 and cols[0] == "MetaData" and "RecordTime" in cols[1]:
308 info["RecordTime"] = cols[2]
309
310 if cols[0] == "DataName":
311 data_started = True
312 continue
313
314 if data_started and cols[0] == "DataValue":
315 if len(cols) < 3:
316 continue
317
318 v = float(cols[1])
319 i = float(cols[2])
320 if not np.isfinite(v) or not np.isfinite(i):
321 continue
322 raw_v.append(v)
323 raw_i.append(i)
324
325 if not raw_v:
326 raise ValueError(f"No valid data found in {infile}")
327
328 V = np.asarray(raw_v, dtype=float)
329 I = np.asarray(raw_i, dtype=float)
330
331 # まず V でソート
332 idx = np.argsort(V)
333 V = V[idx]
334 I = I[idx]
335
336 # 次に xmin/xmax でクリップ
337 mask = np.ones(len(V), dtype=bool)
338 if xmin is not None:
339 mask &= (V >= xmin)
340 if xmax is not None:
341 mask &= (V <= xmax)
342
343 V = V[mask]
344 I = I[mask]
345
346 if len(V) == 0:
347 raise ValueError("No data points remain after xmin/xmax clipping.")
348
349 # 最後に ndataskip を適用
350 if ndataskip > 0:
351 keep = np.zeros(len(V), dtype=bool)
352 keep[::(ndataskip + 1)] = True
353 V = V[keep]
354 I = I[keep]
355
356 print(f"File Name : {info['FileName']}")
357 print(f"Record Time : {info['RecordTime']}")
358 print(f"Points : {len(V)}")
359 print(f"V range : {V.min():.6g} to {V.max():.6g} V")
360
361 return V, I, info
362
363
364def load_param_csv(csv_path):
365 """
366 概要: パラメータ設定CSVファイルからパラメータ値と固定設定を読み込みます。
367
368 詳細説明:
369 CSVファイルには 'varname', 'value', 'optid' の列が含まれると想定されます。
370 'optid' が '0' の場合、そのパラメータは固定されます。
371
372 引数:
373 :param csv_path: 読み込むCSVファイルのパス。
374 :type csv_path: str
375
376 戻り値:
377 :returns: params (読み込まれたパラメータ名と値の辞書), fix_set (固定されるパラメータ名のセット)
378 :rtype: tuple[dict, set]
379 """
380 params = {}
381 fix_set = set()
382
383 if os.path.exists(csv_path):
384 with open(csv_path, "r", encoding="utf-8") as f:
385 reader = csv.DictReader(f)
386 for row in reader:
387 name = row.get("varname", "").strip()
388 if not name:
389 continue
390 try:
391 val = float(row.get("value", ""))
392 except (TypeError, ValueError):
393 continue
394
395 params[name] = val
396 try:
397 if int(row.get("optid", 1)) == 0 and name in PARAM_NAMES:
398 fix_set.add(name)
399 except (TypeError, ValueError):
400 pass
401
402 return params, fix_set
403
404
405def save_param_csv(csv_path, params, fix_set, errors=None, rss=None):
406 """
407 概要: 現在のパラメータ値、固定設定、推定誤差、およびRSSをCSVファイルに保存します。
408
409 詳細説明:
410 PARAM_NAMES にリストされている全てのパラメータと、追加の診断情報（VOC_estなど）が保存されます。
411
412 引数:
413 :param csv_path: 保存するCSVファイルのパス。
414 :type csv_path: str
415 :param params: 保存するパラメータ名と値の辞書。
416 :type params: dict
417 :param fix_set: 固定されたパラメータ名のセット。
418 :type fix_set: set
419 :param errors: 各パラメータの推定誤差を含む辞書。Noneの場合、誤差は0.0として保存。
420 :type errors: dict or None
421 :param rss: 残差平方和 (Residual Sum of Squares)。Noneの場合、保存しない。
422 :type rss: float or None
423
424 戻り値:
425 :returns: なし
426 :rtype: None
427 """
428 if errors is None:
429 errors = {}
430
431 with open(csv_path, "w", encoding="utf-8", newline="") as f:
432 writer = csv.writer(f)
433 writer.writerow(["varname", "value", "optid", "error"])
434 for name in PARAM_NAMES:
435 writer.writerow([
436 name,
437 params[name],
438 0 if name in fix_set else 1,
439 errors.get(name, 0.0),
440 ])
441 for extra in [
442 "VOC_est", "ISC_est", "Vsh_est", "Ish_est",
443 "Vtfe_est", "Itfe_est",
444 ]:
445 if extra in params:
446 writer.writerow([extra, params[extra], 0, 0.0])
447 if rss is not None:
448 writer.writerow(["RSS_logy", rss, 0, 0.0])
449
450def save_iv_to_excel(outfile_xlsx, V, I_input, I_init, I_fit, info,
451 params_init=None, params_final=None, rss=None, errors=None):
452 """
453 概要: IVデータ、初期パラメータ、最終パラメータ、およびサマリー情報をExcelファイルに保存します。
454
455 詳細説明:
456 測定データ、初期モデル電流、フィッティング後のモデル電流が線形および対数絶対値スケールで保存されます。
457 パラメータシートには初期値と最終値、および推定誤差が含まれます。
458
459 引数:
460 :param outfile_xlsx: 出力するExcelファイルのパス。
461 :type outfile_xlsx: str
462 :param V: 電圧データのNumPy配列。
463 :type V: numpy.ndarray
464 :param I_input: 測定された電流データのNumPy配列。Noneの場合、保存しない。
465 :type I_input: numpy.ndarray or None
466 :param I_init: 初期モデルによって計算された電流データのNumPy配列。Noneの場合、保存しない。
467 :type I_init: numpy.ndarray or None
468 :param I_fit: フィッティング後のモデルによって計算された電流データのNumPy配列。Noneの場合、保存しない。
469 :type I_fit: numpy.ndarray or None
470 :param info: ファイル名や記録時間などの情報を含む辞書。
471 :type info: dict
472 :param params_init: 初期パラメータ値を含む辞書。Noneの場合、保存しない。
473 :type params_init: dict or None
474 :param params_final: フィッティング後の最終パラメータ値を含む辞書。Noneの場合、保存しない。
475 :type params_final: dict or None
476 :param rss: 残差平方和 (Residual Sum of Squares)。Noneの場合、保存しない。
477 :type rss: float or None
478 :param errors: 各パラメータの推定誤差を含む辞書。Noneの場合、保存しない。
479 :type errors: dict or None
480
481 戻り値:
482 :returns: なし
483 :rtype: None
484 """
485 print(f"[I/O] Save IV Excel: {outfile_xlsx}")
486
487 wb = Workbook()
488
489 header_fill = PatternFill("solid", fgColor="1F4E78")
490 header_font = Font(color="FFFFFF", bold=True)
491
492 def log_abs_current(x):
493 if x is None:
494 return None
495 return np.log10(np.abs(x) + EPS_I)
496
497 # ------------------------------------------------------------
498 # sheet 1: iv_data
499 # ------------------------------------------------------------
500 ws = wb.active
501 ws.title = "iv_data"
502
503 headers = [
504 "V",
505 "I_input", "I_init", "I_fit",
506 "abs_I_input", "abs_I_init", "abs_I_fit",
507 "log10_abs_I_input", "log10_abs_I_init", "log10_abs_I_fit",
508 ]
509 ws.append(headers)
510 for c in ws[1]:
511 c.fill = header_fill
512 c.font = header_font
513
514 n = len(V)
515 for i in range(n):
516 ii = None if I_input is None else float(I_input[i])
517 i0 = None if I_init is None else float(I_init[i])
518 ifit = None if I_fit is None else float(I_fit[i])
519
520 ws.append([
521 float(V[i]),
522 ii, i0, ifit,
523 None if ii is None else abs(ii),
524 None if i0 is None else abs(i0),
525 None if ifit is None else abs(ifit),
526 None if ii is None else float(log_abs_current(ii)),
527 None if i0 is None else float(log_abs_current(i0)),
528 None if ifit is None else float(log_abs_current(ifit)),
529 ])
530
531 for col in ["A","B","C","D","E","F","G","H","I","J"]:
532 ws.column_dimensions[col].width = 18
533
534 # ------------------------------------------------------------
535 # sheet 2: init_params
536 # ------------------------------------------------------------
537 if params_init is not None:
538 ws2 = wb.create_sheet("init_params")
539 ws2.append(["name", "value"])
540 for c in ws2[1]:
541 c.fill = header_fill
542 c.font = header_font
543 for name in PARAM_NAMES:
544 ws2.append([name, float(params_init.get(name, np.nan))])
545 for extra in ["VOC_est", "ISC_est", "Vsh_est", "Ish_est", "Vtfe_est", "Itfe_est"]:
546 if extra in params_init:
547 ws2.append([extra, float(params_init[extra])])
548
549 # ------------------------------------------------------------
550 # sheet 3: final_params
551 # ------------------------------------------------------------
552 if params_final is not None:
553 ws3 = wb.create_sheet("final_params")
554 ws3.append(["name", "value", "error"])
555 for c in ws3[1]:
556 c.fill = header_fill
557 c.font = header_font
558 if errors is None:
559 errors = {}
560 for name in PARAM_NAMES:
561 ws3.append([
562 name,
563 float(params_final.get(name, np.nan)),
564 float(errors.get(name, np.nan)),
565 ])
566
567 # ------------------------------------------------------------
568 # sheet 4: summary
569 # ------------------------------------------------------------
570 ws4 = wb.create_sheet("summary")
571 ws4.append(["item", "value"])
572 for c in ws4[1]:
573 c.fill = header_fill
574 c.font = header_font
575
576 ws4.append(["FileName", info.get("FileName", "")])
577 ws4.append(["RecordTime", info.get("RecordTime", "")])
578 ws4.append(["Npoints", int(len(V))])
579 if rss is not None:
580 ws4.append(["RSS_logy", float(rss)])
581
582 wb.save(outfile_xlsx)
583 print(f"[I/O] Saved IV Excel: {outfile_xlsx}")
584
585
586# =============================================================================
587# helper math
588# =============================================================================
589
590def solve_root_poly(x, y, target_y=0.0, order=3):
591 """
592 概要: データ点 (x, y) に多項式をフィッティングし、指定された target_y に対応する x の根を近似的に見つけます。
593
594 詳細説明:
595 データ点数が少ない場合は、線形補間または多項式の次数を調整して対処します。
596 複数の実根がある場合は、入力 x の平均に最も近い根を選択します。
597
598 引数:
599 :param x: 独立変数データ点の配列。
600 :type x: numpy.ndarray
601 :param y: 従属変数データ点の配列。
602 :type y: numpy.ndarray
603 :param target_y: 根を見つけたい目標のy値。
604 :type target_y: float
605 :param order: 多項式フィッティングの次数。
606 :type order: int
607
608 戻り値:
609 :returns: target_y に対応する x の推定根。
610 :rtype: float
611 """
612 if len(x) < 2:
613 return float(x[0]) if len(x) else 0.0
614 order = min(order, len(x) - 1)
615 if order < 1:
616 return float(np.interp(target_y, y, x))
617 z = np.polyfit(x, y - target_y, order)
618 roots = np.roots(z)
619 real_roots = roots[np.isreal(roots)].real
620 if len(real_roots) == 0:
621 return float(np.interp(target_y, y, x))
622 return float(real_roots[np.argmin(np.abs(real_roots - np.mean(x)))])
623
624
625def local_poly_value(x, y, x0, npts=7, order=3):
626 """
627 概要: 指定された点 x0 の近傍のデータ点に局所多項式をフィッティングし、x0 におけるy値を推定します。
628
629 詳細説明:
630 npts で指定された数の最も近いデータ点を使用します。
631
632 引数:
633 :param x: 独立変数データ点の配列。
634 :type x: numpy.ndarray
635 :param y: 従属変数データ点の配列。
636 :type y: numpy.ndarray
637 :param x0: y値を推定したい目標のx値。
638 :type x0: float
639 :param npts: 局所フィッティングに使用するデータ点の数。
640 :type npts: int
641 :param order: 多項式フィッティングの次数。
642 :type order: int
643
644 戻り値:
645 :returns: x0 における推定y値。
646 :rtype: float
647 """
648 x = np.asarray(x, dtype=float)
649 y = np.asarray(y, dtype=float)
650 idx = np.argsort(np.abs(x - x0))[:max(2, npts)]
651 xs = x[idx]
652 ys = y[idx]
653 sidx = np.argsort(xs)
654 xs = xs[sidx]
655 ys = ys[sidx]
656 deg = min(order, len(xs) - 1)
657 if deg < 1: # データ点が1点のみの場合、その点のy値を返す
658 if len(xs) == 1:
659 return float(ys[0])
660 else: # データ点が2点の場合、線形補間
661 return float(np.interp(x0, xs, ys))
662
663 coef = np.polyfit(xs, ys, deg)
664 return float(np.polyval(coef, x0))
665
666
667def smooth_polyfit(y, window_points=5, poly_order=3):
668 """
669 概要: Savitzky-Golayフィルターに似た局所多項式フィッティングにより、配列 y を平滑化します。
670
671 詳細説明:
672 各データ点に対して、その点と近傍の window_points 数（奇数）のデータ点に
673 poly_order 次の多項式をフィッティングし、中央の点の値を推定します。
674
675 引数:
676 :param y: 平滑化するデータ点の配列。
677 :type y: numpy.ndarray
678 :param window_points: 各局所フィッティングに使用するウィンドウの点の数（奇数であるべき）。
679 :type window_points: int
680 :param poly_order: 各局所フィッティングで使用する多項式の次数。
681 :type poly_order: int
682
683 戻り値:
684 :returns: 平滑化されたデータ点の配列。
685 :rtype: numpy.ndarray
686 """
687 y = np.asarray(y, dtype=float)
688 n = len(y)
689
690 if n < 3:
691 return y.copy()
692
693 if window_points < 3:
694 window_points = 3
695 if window_points % 2 == 0:
696 window_points += 1
697 if window_points > n:
698 window_points = n if (n % 2 == 1) else n - 1
699 if window_points < 3: # ここに到達したら、平滑化不能なのでコピーを返す
700 return y.copy()
701
702 if poly_order >= window_points:
703 poly_order = window_points - 1
704 if poly_order < 1:
705 poly_order = 1
706
707 half = window_points // 2
708 ys = np.empty(n, dtype=float)
709
710 for i in range(n):
711 i0 = max(0, i - half)
712 i1 = min(n, i + half + 1)
713 while (i1 - i0) < window_points:
714 if i0 > 0:
715 i0 -= 1
716 elif i1 < n:
717 i1 += 1
718 else:
719 break
720
721 idx = np.arange(i0, i1, dtype=float)
722 yy = y[i0:i1]
723 xloc = idx - i # 中心を0とする相対座標
724
725 deg = min(poly_order, len(yy) - 1)
726 if deg < 1:
727 ys[i] = y[i]
728 continue
729
730 coeff = np.polyfit(xloc, yy, deg)
731 ys[i] = np.polyval(coeff, 0.0)
732
733 return ys
734
735
736def sigmoid_blend(x, dV):
737 """
738 概要: シグモイド関数を用いて、2つのモデル間のブレンド係数を計算します。
739
740 詳細説明:
741 x が0.0の場合に0.5、x が正で十分に大きい場合に1.0、
742 x が負で十分に小さい場合に0.0に近づく値を返します。
743 dV はブレンドの遷移幅を決定します。dV が非常に小さい場合、ステップ関数に近づきます。
744
745 引数:
746 :param x: ブレンドの基準となる電圧または汎用的な値。
747 :type x: numpy.ndarray or float
748 :param dV: シグモイド遷移の幅を制御するパラメータ。
749 :type dV: float
750
751 戻り値:
752 :returns: x に応じたブレンド係数。0から1の間の値を取ります。
753 :rtype: numpy.ndarray or float
754 """
755 x = np.asarray(x, dtype=float)
756 if dV < 1.0e-6:
757 return (x >= 0.0).astype(float)
758 z = np.clip(-x / dV, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
759 return 1.0 / (1.0 + np.exp(z))
760
761
762# =============================================================================
763# current components
764# =============================================================================
765
766def j_diode(Vd, I0, ndiode, temperature):
767 """
768 概要: 理想ダイオードモデル（Shockley方程式）の電流成分を計算します。
769
770 引数:
771 :param Vd: ダイオード接合にかかる電圧。
772 :type Vd: numpy.ndarray or float
773 :param I0: 逆飽和電流。
774 :type I0: float
775 :param ndiode: ダイオード理想因子（エミッタ係数）。
776 :type ndiode: float
777 :param temperature: 温度 (K)。
778 :type temperature: float
779
780 戻り値:
781 :returns: ダイオード電流。
782 :rtype: numpy.ndarray or float
783 """
784 vt = ndiode * KB * temperature / E_CHARGE
785 arg = np.clip(Vd / vt, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
786 return I0 * (np.exp(arg) - 1.0)
787
788
789def j_tfe_forward(Vd, A_tfe, E00):
790 """
791 概要: 熱アシストトンネルフィールドエミッション (TFE) 電流のフォワードバイアス成分を計算します。
792
793 詳細説明:
794 TFEは通常、低温や欠陥の多い材料で支配的になるメカニズムです。
795 この関数は、Vd > 0 の領域でTFEによる電流を計算します。
796
797 引数:
798 :param Vd: ダイオード接合にかかる電圧。
799 :type Vd: numpy.ndarray or float
800 :param A_tfe: TFE電流の振幅係数。
801 :type A_tfe: float
802 :param E00: TFEの特性エネルギーまたは係数。
803 :type E00: float
804
805 戻り値:
806 :returns: TFEフォワード電流。
807 :rtype: numpy.ndarray or float
808 """
809 e00 = max(E00, 1.0e-12) # ゼロ除算防止
810 arg = np.clip(Vd / e00, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
811 return A_tfe * (np.exp(arg) - 1.0)
812
813
814def j_tfe_reverse(Vd, A_tfe, E00):
815 """
816 概要: 熱アシストトンネルフィールドエミッション (TFE) 電流のリバースバイアス成分を計算します。
817
818 詳細説明:
819 この関数は、Vd < 0 の領域でTFEによる電流を計算します。
820
821 引数:
822 :param Vd: ダイオード接合にかかる電圧。
823 :type Vd: numpy.ndarray or float
824 :param A_tfe: TFE電流の振幅係数。
825 :type A_tfe: float
826 :param E00: TFEの特性エネルギーまたは係数。
827 :type E00: float
828
829 戻り値:
830 :returns: TFEリバース電流。
831 :rtype: numpy.ndarray or float
832 """
833 e00 = max(E00, 1.0e-12) # ゼロ除算防止
834 arg = np.clip(-Vd / e00, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
835 return -A_tfe * (np.exp(arg) - 1.0)
836
837
838def j_fn_forward(Vd, A_fn, B_fn):
839 """
840 概要: Fowler-Nordheimトンネル電流 (FN) のフォワードバイアス成分を計算します。
841
842 詳細説明:
843 このメカニズムは、高注入レベルや特定の欠陥メカニズムによって引き起こされる電流をモデル化するために使用できます。
844 Vd > 0 の領域で定義されます。
845
846 引数:
847 :param Vd: ダイオード接合にかかる電圧。
848 :type Vd: numpy.ndarray or float
849 :param A_fn: FN電流の振幅係数。
850 :type A_fn: float
851 :param B_fn: FN電流の指数関数的電圧依存性を制御する係数。
852 :type B_fn: float
853
854 戻り値:
855 :returns: FNフォワード電流。
856 :rtype: numpy.ndarray or float
857 """
858 if np.ndim(Vd) == 0:
859 if Vd <= 1.0e-12:
860 return 0.0
861 arg = np.clip(-B_fn / Vd, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
862 return A_fn * (Vd ** 2) * np.exp(arg)
863
864 arr = np.asarray(Vd, dtype=float)
865 out = np.zeros_like(arr)
866 mask = arr > 1.0e-12
867 arg = np.clip(-B_fn / arr[mask], -700.0, 700.0) # オーバーフロー防止のためのクリッピング
868 out[mask] = A_fn * (arr[mask] ** 2) * np.exp(arg)
869 return out
870
871
872def j_fn_reverse(Vd, A_fn, B_fn):
873 """
874 概要: Fowler-Nordheimトンネル電流 (FN) のリバースバイアス成分を計算します。
875
876 詳細説明:
877 このメカニズムは、高注入レベルや特定の欠陥メカニズムによって引き起こされる電流をモデル化するために使用できます。
878 Vd < 0 の領域で定義されます。
879
880 引数:
881 :param Vd: ダイオード接合にかかる電圧。
882 :type Vd: numpy.ndarray or float
883 :param A_fn: FN電流の振幅係数。
884 :type A_fn: float
885 :param B_fn: FN電流の指数関数的電圧依存性を制御する係数。
886 :type B_fn: float
887
888 戻り値:
889 :returns: FNリバース電流。
890 :rtype: numpy.ndarray or float
891 """
892 if np.ndim(Vd) == 0:
893 if Vd >= -1.0e-12:
894 return 0.0
895 vabs = -Vd
896 arg = np.clip(-B_fn / vabs, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
897 return -A_fn * (vabs ** 2) * np.exp(arg)
898
899 arr = np.asarray(Vd, dtype=float)
900 out = np.zeros_like(arr)
901 mask = arr < -1.0e-12
902 vabs = -arr[mask]
903 arg = np.clip(-B_fn / vabs, -700.0, 700.0) # オーバーフロー防止のためのクリッピング
904 out[mask] = -A_fn * (vabs ** 2) * np.exp(arg)
905 return out
906
907
908def j_sclc_trap_transition(Vd, K_sclc, m_sclc, Vtfl, dV=1e-3):
909 """
910 概要: トラップ制限空間電荷制限電流（SCLC）とトラップフリーSCLCの間の遷移をモデル化する電流成分を計算します。
911
912 詳細説明:
913 このモデルは、電圧に応じてトラップが満たされることで、電流が異なる電圧依存性を示す領域を表現します。
914 具体的には、トラップ制限領域では Vd^m_sclc に、トラップフリー領域では Vd^2 に従います。
915 シグモイド関数を用いて両者間の滑らかな遷移を実現します。
916
917 引数:
918 :param Vd: ダイオード接合にかかる電圧。
919 :type Vd: numpy.ndarray or float
920 :param K_sclc: SCLCの振幅係数（トラップ制限領域の係数）。
921 :type K_sclc: float
922 :param m_sclc: トラップ制限領域でのSCLCの電圧依存性指数。
923 :type m_sclc: float
924 :param Vtfl: トラップが満たされる閾値電圧。
925 :type Vtfl: float
926 :param dV: 遷移の幅を制御する電圧（シグモイド関数の傾き）。
927 :type dV: float
928
929 戻り値:
930 :returns: トラップ遷移を考慮したSCLC電流。
931 :rtype: numpy.ndarray or float
932 """
933 vabs = np.abs(Vd)
934 K_tf = K_sclc * (max(Vtfl, 1.0e-30) ** (m_sclc - 2.0))
935
936 if dV < 1.0e-12:
937 w = (vabs >= Vtfl).astype(float) if np.ndim(vabs) else float(vabs >= Vtfl)
938 else:
939 z = np.clip(-(vabs - Vtfl) / dV, -700.0, 700.0)
940 w = 1.0 / (1.0 + np.exp(z))
941
942 i_low = K_sclc * (vabs ** m_sclc) # trap-limited
943 i_high = K_tf * (vabs ** 2.0) # trap-free
944 imag = (1.0 - w) * i_low + w * i_high
945 return np.sign(Vd) * imag
946
947def j_sclc(Vd, K_sclc):
948 """
949 概要: 空間電荷制限電流 (SCLC) を計算します。
950
951 詳細説明:
952 SCLCは、高抵抗材料におけるキャリア注入が空間電荷によって制限されるときに発生します。
953 この電流は電圧の2乗に比例し、電圧の符号によって方向が決まります。
954
955 引数:
956 :param Vd: ダイオード接合にかかる電圧。
957 :type Vd: numpy.ndarray or float
958 :param K_sclc: SCLCの振幅係数。
959 :type K_sclc: float
960
961 戻り値:
962 :returns: SCLC電流。
963 :rtype: numpy.ndarray or float
964 """
965 return np.sign(Vd) * K_sclc * (np.abs(Vd) ** 2)
966
967
968def parse_mech_mode(mode_str):
969 """
970 概要: 輸送メカニズムのモード文字列を解析し、使用するメカニズムのセットを返します。
971
972 詳細説明:
973 モード文字列はカンマ (,) またはプラス (+) で区切られたメカニズム名（例: diode+tfe）を
974 含みます。無効なメカニズム名が指定された場合はValueErrorを発生させます。
975
976 引数:
977 :param mode_str: 解析するモード文字列（例: diode, diode+tfe, none）。
978 :type mode_str: str or None
979
980 戻り値:
981 :returns: 使用するメカニズム名のセット。
982 :rtype: set
983
984 例外:
985 :raises ValueError: 未知の輸送メカニズム名が指定された場合。
986 """
987 if mode_str is None:
988 return set()
989 txt = mode_str.lower().replace(" ", "")
990 if txt in ("", "none", "off"):
991 return set()
992 parts = [p for p in txt.replace(",", "+").split("+") if p]
993 norm = set()
994 for p in parts:
995 if p in ("diode", "tfe", "fn", "sclc"):
996 norm.add(p)
997 else:
998 raise ValueError(f"Unknown transport mechanism: {p}")
999 return norm
1000
1001
1002def used_params_from_modes(forwardIV, reverseIV):
1003 """
1004 概要: 指定されたフォワードおよびリバースIVモデルで使用されるパラメータ名を決定します。
1005
1006 詳細説明:
1007 共通のパラメータ (Rs, Rsh, IPV) に加えて、
1008 有効な輸送メカニズム (diode, tfe, fn, sclc) に関連するパラメータを特定します。
1009
1010 引数:
1011 :param forwardIV: フォワードバイアスモデルのメカニズムを示す文字列（例: diode, diode+tfe）。
1012 :type forwardIV: str
1013 :param reverseIV: リバースバイアスモデルのメカニズムを示す文字列（例: diode+tfe, diode+fn）。
1014 :type reverseIV: str
1015
1016 戻り値:
1017 :returns: 使用されるパラメータ名のセット。
1018 :rtype: set
1019 """
1020 used = {"Rs", "Rsh", "IPV"}
1021 all_modes = parse_mech_mode(forwardIV) | parse_mech_mode(reverseIV)
1022 if "diode" in all_modes:
1023 used |= {"I0", "ndiode"}
1024 if "tfe" in all_modes:
1025 used |= {"A_tfe", "E00"}
1026 if "fn" in all_modes:
1027 used |= {"A_fn", "B_fn"}
1028 if "sclc" in all_modes:
1029 used |= {"K_sclc"}
1030 return used
1031
1032
1033def junction_current_forward(Vd, p, mode_set, temperature):
1034 """
1035 概要: フォワードバイアス条件における総接合電流を計算します。
1036
1037 詳細説明:
1038 mode_set で指定されたメカニズムに基づいて、対応する電流成分（ダイオード、TFE、FN、SCLC）を合計します。
1039
1040 引数:
1041 :param Vd: ダイオード接合にかかる電圧。
1042 :type Vd: numpy.ndarray or float
1043 :param p: パラメータ値を含む辞書。
1044 :type p: dict
1045 :param mode_set: フォワードバイアスで使用するメカニズムのセット。
1046 :type mode_set: set
1047 :param temperature: 温度 (K)。
1048 :type temperature: float
1049
1050 戻り値:
1051 :returns: 総接合電流。
1052 :rtype: numpy.ndarray or float
1053 """
1054 j = 0.0
1055 if "diode" in mode_set:
1056 j += j_diode(Vd, p["I0"], p["ndiode"], temperature)
1057 if "tfe" in mode_set:
1058 j += j_tfe_forward(Vd, p["A_tfe"], p["E00"])
1059 if "fn" in mode_set:
1060 j += j_fn_forward(Vd, p["A_fn"], p["B_fn"])
1061 if "sclc" in mode_set:
1062 j += j_sclc(Vd, p["K_sclc"])
1063 return j
1064
1065
1066def junction_current_reverse(Vd, p, mode_set, temperature):
1067 """
1068 概要: リバースバイアス条件における総接合電流を計算します。
1069
1070 詳細説明:
1071 mode_set で指定されたメカニズムに基づいて、対応する電流成分（ダイオード、TFE、FN、SCLC）を合計します。
1072
1073 引数:
1074 :param Vd: ダイオード接合にかかる電圧。
1075 :type Vd: numpy.ndarray or float
1076 :param p: パラメータ値を含む辞書。
1077 :type p: dict
1078 :param mode_set: リバースバイアスで使用するメカニズムのセット。
1079 :type mode_set: set
1080 :param temperature: 温度 (K)。
1081 :type temperature: float
1082
1083 戻り値:
1084 :returns: 総接合電流。
1085 :rtype: numpy.ndarray or float
1086 """
1087 j = 0.0
1088 if "diode" in mode_set:
1089 j += j_diode(Vd, p["I0"], p["ndiode"], temperature)
1090 if "tfe" in mode_set:
1091 j += j_tfe_reverse(Vd, p["A_tfe"], p["E00"])
1092 if "fn" in mode_set:
1093 j += j_fn_reverse(Vd, p["A_fn"], p["B_fn"])
1094 if "sclc" in mode_set:
1095 j += j_sclc(Vd, p["K_sclc"])
1096 return j
1097
1098
1099# =============================================================================
1100# model
1101# =============================================================================
1102
1103def same_sign(a, b):
1104 """
1105 概要: 2つの数値が同じ符号を持つかどうかを判定します。
1106
1107 詳細説明:
1108 両方が正または両方が負の場合にTrueを返します。
1109 片方または両方がゼロの場合の動作は、厳密には考慮されていません。
1110
1111 引数:
1112 :param a: 1つ目の数値。
1113 :type a: float
1114 :param b: 2つ目の数値。
1115 :type b: float
1116
1117 戻り値:
1118 :returns: a と b が同じ符号を持つ場合はTrue、そうでない場合はFalse。
1119 :rtype: bool
1120 """
1121 return (a > 0.0 and b > 0.0) or (a < 0.0 and b < 0.0)
1122
1123def model(V, p, forwardIV="diode", reverseIV="diode+tfe", dV=0.05, temperature=300.0):
1124 """
1125 概要: 拡張一ダイオードモデル（SDM）に基づいて電流-電圧（IV）カーブを計算します。
1126
1127 詳細説明:
1128 与えられた外部電圧 V、パラメータ p、および選択された輸送メカニズム (forwardIV, reverseIV) を使用して
1129 対応する電流 I を計算します。接合電圧 Vd を未知数とする根探しアルゴリズム（Brentのメソッド）を
1130 用いて安定的にモデル電流を求めます。
1131 フォワードとリバースで異なるモデルが指定された場合、シグモイド関数 (sigmoid_blend) を使って
1132 Vd が0V付近でスムーズにブレンドされます。
1133
1134 引数:
1135 :param V: 外部印加電圧の配列。
1136 :type V: numpy.ndarray
1137 :param p: モデルパラメータを含む辞書。主要なキーは I0, ndiode, IPV, Rs, Rsh, A_tfe, E00, A_fn, B_fn, K_sclc。
1138 :type p: dict
1139 :param forwardIV: フォワードバイアス領域で使用する輸送メカニズムを示す文字列。例: diode, diode+tfe+sclc。
1140 :type forwardIV: str
1141 :param reverseIV: リバースバイアス領域で使用する輸送メカニズムを示す文字列。例: diode+tfe, diode+fn。
1142 :type reverseIV: str
1143 :param dV: フォワードモデルとリバースモデルをブレンドする際の電圧幅。Vd=0V付近での遷移の鋭さを制御します。
1144 :type dV: float
1145 :param temperature: シミュレーションの温度 (K)。
1146 :type temperature: float
1147
1148 戻り値:
1149 :returns: 各外部電圧 V に対応する総電流の配列。
1150 :rtype: numpy.ndarray
1151 """
1152 fset = parse_mech_mode(forwardIV)
1153 rset = parse_mech_mode(reverseIV)
1154 same_model = (fset == rset)
1155
1156 V = np.atleast_1d(np.asarray(V, dtype=float))
1157 Itotal = np.empty_like(V)
1158
1159 Rs = max(float(p["Rs"]), EPS_R)
1160 Rsh = max(float(p["Rsh"]), EPS_R)
1161 IPV = float(p["IPV"])
1162
1163 def residual_vd(vd_trial, v_target):
1164 """
1165 概要: 接合電圧 vd_trial の残差関数を計算します。
1166 この関数は brentq によって vd_trial の根を見つけるために使用されます。
1167 """
1168 if same_model:
1169 ij = junction_current_forward(vd_trial, p, fset, temperature)
1170 else:
1171 w = float(sigmoid_blend(vd_trial, dV))
1172 ij_f = junction_current_forward(vd_trial, p, fset, temperature)
1173 ij_r = junction_current_reverse(vd_trial, p, rset, temperature)
1174 ij = w * ij_f + (1.0 - w) * ij_r
1175
1176 ibranch = ij + vd_trial / Rsh
1177 v_calc = vd_trial + Rs * ibranch
1178 return v_calc - v_target
1179
1180 for i, v_val in enumerate(V):
1181 try:
1182 # 根探しの初期探索範囲を設定。VdはVよりRs*I落ちるため、広めに設定。
1183 width = max(5.0, abs(v_val) + Rs * (abs(IPV) + abs(v_val) / Rsh + 1.0))
1184 lo = v_val - width
1185 hi = v_val + width
1186
1187 # BrentQ のための探索範囲 [lo, hi] を確保。f(lo) と f(hi) が異なる符号を持つように調整。
1188 flo = residual_vd(lo, v_val)
1189 fhi = residual_vd(hi, v_val)
1190 ntry = 0
1191 while same_sign(flo, fhi) and ntry < 25:
1192 width *= 2.0
1193 lo = v_val - width
1194 hi = v_val + width
1195 flo = residual_vd(lo, v_val)
1196 fhi = residual_vd(hi, v_val)
1197 ntry += 1
1198
1199 if same_sign(flo, fhi): # 符号が異ならない場合、根が見つからなかったとみなし、近似値を返す
1200 vd_root = v_val
1201 else:
1202 vd_root = brentq(
1203 residual_vd, lo, hi,
1204 args=(v_val,),
1205 maxiter=300, xtol=1.0e-12, rtol=1.0e-10
1206 )
1207
1208 # 根が見つかったVdから電流を再計算
1209 if same_model:
1210 ij = junction_current_forward(vd_root, p, fset, temperature)
1211 else:
1212 w = float(sigmoid_blend(vd_root, dV))
1213 ij_f = junction_current_forward(vd_root, p, fset, temperature)
1214 ij_r = junction_current_reverse(vd_root, p, rset, temperature)
1215 ij = w * ij_f + (1.0 - w) * ij_r
1216
1217 ibranch = ij + vd_root / Rsh
1218 Itotal[i] = ibranch - IPV
1219 except Exception:
1220 # 根探しに失敗した場合のフォールバック（例：V = 0 で Rs + Rsh がゼロに近い場合など）
1221 Itotal[i] = v_val / max(Rs + Rsh, EPS_R) - IPV
1222
1223 return Itotal
1224
1225
1226# =============================================================================
1227# initial parameter estimation
1228# =============================================================================
1229
1230def estimate_voc_local_quadratic(V, I, nlsq_points=5):
1231 """
1232 概要: IVデータから開回路電圧（VOC）を局所二次多項式フィッティングで推定します。
1233
1234 詳細説明:
1235 電流がゼロに最も近い点の近傍データを使用して、二次多項式をフィッティングし、
1236 その多項式の根からVOCを求めます。データ点が少ない場合は、線形補間を使用します。
1237
1238 引数:
1239 :param V: 電圧データの配列。
1240 :type V: numpy.ndarray
1241 :param I: 電流データの配列。
1242 :type I: numpy.ndarray
1243 :param nlsq_points: 局所フィッティングに使用するデータ点の数。
1244 :type nlsq_points: int
1245
1246 戻り値:
1247 :returns: 推定された開回路電圧VOC。
1248 :rtype: float
1249 """
1250 V = np.asarray(V, dtype=float)
1251 I = np.asarray(I, dtype=float)
1252
1253 if len(V) < 3:
1254 return float(V[np.argmin(np.abs(I))])
1255
1256 nlsq_points = max(3, int(nlsq_points))
1257 nlsq_points = min(nlsq_points, len(V))
1258
1259 ic = int(np.argmin(np.abs(I))) # abs(I) 最小点
1260 half = nlsq_points // 2
1261
1262 i0 = max(0, ic - half)
1263 i1 = min(len(V), i0 + nlsq_points)
1264 i0 = max(0, i1 - nlsq_points)
1265
1266 x = V[i0:i1]
1267 y = I[i0:i1]
1268
1269 deg = min(2, len(x) - 1)
1270 coef = np.polyfit(x, y, deg)
1271
1272 if deg == 1:
1273 a, b = coef
1274 if abs(a) < 1.0e-30:
1275 return float(V[ic])
1276 return float(-b / a)
1277
1278 a, b, c = coef
1279 if abs(a) < 1.0e-30:
1280 if abs(b) < 1.0e-30:
1281 return float(V[ic])
1282 return float(-c / b)
1283
1284 roots = np.roots([a, b, c])
1285 real_roots = roots[np.isreal(roots)].real
1286
1287 if len(real_roots) == 0:
1288 return float(V[ic])
1289
1290 # 近傍窓の中心に最も近い根を採用
1291 x_center = V[ic]
1292 return float(real_roots[np.argmin(np.abs(real_roots - x_center))])
1293
1294def estimate_initial_params(V, I_meas, temperature):
1295 """
1296 概要: IVデータからSDM（拡張一ダイオードモデル）の初期パラメータを推定します。
1297
1298 詳細説明:
1299 データ平滑化、dI/dV解析、多項式根探し、負電圧領域でのTFE近似を用いて、
1300 最適化に適した頑健な初期値を導出します。
1301
1302 引数:
1303 :param V: 測定された電圧データ点の配列。
1304 :type V: numpy.ndarray
1305 :param I_meas: 測定された電流データ点の配列。
1306 :type I_meas: numpy.ndarray
1307 :param temperature: 測定時の温度 (K)。
1308 :type temperature: float
1309
1310 戻り値:
1311 :returns: 推定された初期パラメータを含む辞書。
1312 :rtype: dict
1313 """
1314 print("[Process] Estimating initial diode/TFE parameters from data...")
1315
1316 # 重複するV値を処理
1317 uV, idx = np.unique(V, return_index=True)
1318 uI = I_meas[idx]
1319
1320 # データ平滑化と微分
1321 y_sm = smooth_polyfit(uI, window_points=5, poly_order=3)
1322 dydv = np.gradient(y_sm, uV)
1323
1324 # Rsの推定: V>0 で dI/dV が最大となる点から
1325 pos_mask = uV > 0.0
1326 if np.any(pos_mask):
1327 dpos = dydv[pos_mask]
1328 xpos = uV[pos_mask]
1329 idx_pos_local = int(np.argmax(dpos))
1330 V_rs = float(xpos[idx_pos_local])
1331 slope_pos = float(dpos[idx_pos_local])
1332 else: # V>0 のデータがない場合、全体の最大勾配
1333 V_rs = float(uV[-1])
1334 slope_pos = float(np.max(dydv))
1335 Rs_est = 1.0 / abs(slope_pos)
1336
1337 # Rshの推定: V<0 で dI/dV が最小（最も負）となる点から
1338 neg_mask = uV < 0.0
1339 if np.any(neg_mask):
1340 dneg = dydv[neg_mask]
1341 xneg = uV[neg_mask]
1342 idx_neg_local = int(np.argmin(dneg))
1343 V_sh = float(xneg[idx_neg_local])
1344 slope_neg = float(dneg[idx_neg_local])
1345 else: # V<0 のデータがない場合、全体の最小勾配
1346 V_sh = float(uV[0])
1347 slope_neg = float(np.min(dydv))
1348 Rsh_est = 1.0 / max(abs(slope_neg), 1.0e-12)
1349
1350 # ISC, IPV, VOC の推定
1351 ISC_est = local_poly_value(uV, y_sm, 0.0, npts=7, order=3)
1352 IPV_est = max(-ISC_est, EPS_I)
1353 VOC_est = solve_root_poly(uV, y_sm, target_y=0.0, order=3)
1354 VOC_est = estimate_voc_local_quadratic(uV, y_sm, nlsq_points=5)
1355
1356 # Ish, I0 の推定
1357 Ish_est = local_poly_value(uV, y_sm, V_sh, npts=7, order=3)
1358 I0_est = max(ISC_est - Ish_est, 1.0e-12) # I0は通常正の値
1359
1360 # TFE 初期値推定 (V<0 の十分負側で ln(abs(I))-V を一次近似)
1361 Vtfe_est = np.nan
1362 Itfe_est = np.nan
1363 A_tfe_est = 1.0e-30
1364 E00_est = 1.0e-3
1365
1366 neg_idx = np.where(uV < 0.0)[0]
1367 if len(neg_idx) >= 5: # 負電圧領域に十分なデータ点がある場合
1368 nneg = len(neg_idx)
1369 ntail = max(5, int(np.ceil(0.4 * nneg))) # 最も負側の40%または5点以上を使用
1370 tail_idx = neg_idx[:ntail] # 最も負側の領域
1371
1372 V_tail = uV[tail_idx]
1373 I_tail = np.abs(y_sm[tail_idx])
1374
1375 mask_good = np.isfinite(V_tail) & np.isfinite(I_tail) & (I_tail > EPS_I)
1376 V_tail = V_tail[mask_good]
1377 I_tail = I_tail[mask_good]
1378
1379 if len(V_tail) >= 2: # 少なくとも2点あれば直線フィッティング可能
1380 lnI = np.log(I_tail)
1381 slope_tfe, intercept_tfe = np.polyfit(V_tail, lnI, 1)
1382
1383 if abs(slope_tfe) > 1.0e-12:
1384 E00_est = 1.0 / abs(slope_tfe)
1385 else:
1386 E00_est = 1.0e-3 # デフォルト値
1387
1388 irep = int(np.argmin(V_tail)) # 最も負の電圧点
1389 Vtfe_est = float(V_tail[irep])
1390 Itfe_est = float(I_tail[irep])
1391
1392 # ln|I| = ln(A_tfe) + (-V)/E00 for reverse-TFE-like growth
1393 # |I| = A_tfe * exp(-V/E00) より
1394 A_tfe_est = Itfe_est * np.exp(Vtfe_est / max(E00_est, 1.0e-12))
1395
1396 A_tfe_est = max(float(A_tfe_est), 1.0e-30)
1397 E00_est = max(float(E00_est), 1.0e-6)
1398
1399 # 推定されたパラメータとデフォルト値を結合
1400 return {
1401 "I0": I0_est,
1402 "ndiode": 1.5,
1403 "IPV": IPV_est,
1404 "Rs": max(Rs_est, 1.0e-6), # Rsはゼロにならないようにする
1405 "Rsh": max(Rsh_est, 1.0e3), # Rshは十分に大きい値から始める
1406 "A_tfe": A_tfe_est,
1407 "E00": E00_est,
1408 "A_fn": 1.0e-30, # FN/SCLCは初期値推定でカバーされないため極小値
1409 "B_fn": 10.0,
1410 "K_sclc": 1.0e-30,
1411 "VOC_est": VOC_est,
1412 "ISC_est": ISC_est,
1413 "Vsh_est": V_sh,
1414 "Ish_est": Ish_est,
1415 "Vtfe_est": Vtfe_est,
1416 "Itfe_est": Itfe_est,
1417 "dIdV_pos_rep": slope_pos,
1418 "dIdV_neg_rep": slope_neg,
1419 }
1420
1421
1422# =============================================================================
1423# fit support
1424# =============================================================================
1425
1426def pack_free_params(params, fix_set):
1427 """
1428 概要: 辞書形式の全パラメータから、最適化対象の自由パラメータのみをNumPy配列にパックします。
1429
1430 詳細説明:
1431 fix_set に含まれるパラメータは除外されます。
1432 LOG_PARAMS に含まれるパラメータは log10 変換されて配列に格納されます。
1433
1434 引数:
1435 :param params: 全パラメータ（固定および自由）を含む辞書。
1436 :type params: dict
1437 :param fix_set: 固定されるパラメータ名のセット。
1438 :type fix_set: set
1439
1440 戻り値:
1441 :returns: 最適化に使用される自由パラメータのNumPy配列。
1442 :rtype: numpy.ndarray
1443 """
1444 p0 = []
1445 for name in PARAM_NAMES:
1446 if name in fix_set:
1447 continue
1448 val = float(params[name])
1449 if name in LOG_PARAMS:
1450 p0.append(np.log10(max(val, 1.0e-30))) # ゼロや負の値を避ける
1451 else:
1452 p0.append(val)
1453 return np.asarray(p0, dtype=float)
1454
1455
1456def unpack_free_params(p_free, base_params, fix_set):
1457 """
1458 概要: 自由パラメータのNumPy配列を、固定パラメータと組み合わせて辞書形式の全パラメータにアンパックします。
1459
1460 詳細説明:
1461 base_params は固定パラメータの値を決定するために使用されます。
1462 LOG_PARAMS に含まれるパラメータは 10** 変換されて辞書に格納されます。
1463
1464 引数:
1465 :param p_free: 最適化中の自由パラメータのNumPy配列。
1466 :type p_free: numpy.ndarray
1467 :param base_params: 固定パラメータの基準値を含む辞書。
1468 :type base_params: dict
1469 :param fix_set: 固定されたパラメータ名のセット。
1470 :type fix_set: set
1471
1472 戻り値:
1473 :returns: 更新された全パラメータを含む辞書。
1474 :rtype: dict
1475 """
1476 curr = dict(base_params) # base_params をコピーして、固定されていないパラメータを更新
1477 idx = 0
1478 for name in PARAM_NAMES:
1479 if name in fix_set:
1480 continue
1481 val = float(p_free[idx])
1482 if name in LOG_PARAMS:
1483 curr[name] = 10.0 ** val
1484 else:
1485 curr[name] = val
1486 idx += 1
1487 return curr
1488
1489
1490
1491def signed_log10_current(I, eps=EPS_I):
1492 """
1493 概要: 符号を保持したまま電流を圧縮する補助関数。
1494
1495 詳細説明:
1496 sign(I) * log10(abs(I)+eps) を返します。
1497 電流値がゼロに近い場合に計算が不安定になることを防ぐため、abs(I) に eps を加えます。
1498
1499 引数:
1500 :param I: 電流値の配列または単一値。
1501 :type I: numpy.ndarray or float
1502 :param eps: 絶対値計算時に加算される微小量。ゼロ除算や対数関数の引数が負になるのを防ぐ。
1503 :type eps: float
1504
1505 戻り値:
1506 :returns: 符号を保持したログスケールの電流値。
1507 :rtype: numpy.ndarray or float
1508 """
1509 I = np.asarray(I, dtype=float)
1510 return np.sign(I) * np.log10(np.abs(I) + eps)
1511
1512
1513def objective(p_free, V, I_meas, args, base_params, fix_set):
1514 """
1515 概要: 最適化のための目的関数（残差平方和）を計算します。
1516
1517 詳細説明:
1518 測定電流とモデル電流の符号を保持した圧縮表現
1519 sign(I) * log10(abs(I)+EPS_I) の残差平方和を使用します。
1520
1521 引数:
1522 :param p_free: 最適化中の自由パラメータのNumPy配列。
1523 :type p_free: numpy.ndarray
1524 :param V: 測定された電圧データ点の配列。
1525 :type V: numpy.ndarray
1526 :param I_meas: 測定された電流データ点の配列。
1527 :type I_meas: numpy.ndarray
1528 :param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト（例: forwardIV, reverseIV）。
1529 :type args: argparse.Namespace
1530 :param base_params: 固定パラメータの基準値を含む辞書。
1531 :type base_params: dict
1532 :param fix_set: 固定されたパラメータ名のセット。
1533 :type fix_set: set
1534
1535 戻り値:
1536 :returns: 目的関数の値（残差平方和）。
1537 :rtype: float
1538 """
1539 curr = unpack_free_params(p_free, base_params, fix_set)
1540 I_calc = model(
1541 V, curr,
1542 forwardIV=args.forwardIV,
1543 reverseIV=args.reverseIV,
1544 dV=args.dv_blend,
1545 temperature=args.temperature,
1546 )
1547 # 符号付き log 圧縮スケールでの残差平方和を計算。
1548 y_meas = signed_log10_current(I_meas, EPS_I)
1549 y_calc = signed_log10_current(I_calc, EPS_I)
1550 return float(np.sum((y_meas - y_calc) ** 2))
1551
1552
1553# =============================================================================
1554# error estimation
1555# =============================================================================
1556
1557def get_jacobian(p_free, V, temperature, args, base_params, fix_set):
1558 """
1559 概要: モデルのヤコビ行列を数値的に計算します。
1560
1561 詳細説明:
1562 各自由パラメータを微小量 eps だけ摂動させ、log10(abs(I)) の変化を計算して
1563 偏微分（ヤコビ行列の要素）を近似します。
1564
1565 引数:
1566 :param p_free: 最適化された自由パラメータのNumPy配列。
1567 :type p_free: numpy.ndarray
1568 :param V: 測定された電圧データ点の配列。
1569 :type V: numpy.ndarray
1570 :param temperature: シミュレーションの温度 (K)。
1571 :type temperature: float
1572 :param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト。
1573 :type args: argparse.Namespace
1574 :param base_params: 固定パラメータの基準値を含む辞書。
1575 :type base_params: dict
1576 :param fix_set: 固定されたパラメータ名のセット。
1577 :type fix_set: set
1578
1579 戻り値:
1580 :returns: 各データ点に対する各自由パラメータの偏微分からなるヤコビ行列。
1581 :rtype: numpy.ndarray
1582 """
1583 eps = 1.0e-5 # 微小摂動量
1584 n_params = len(p_free)
1585 n_points = len(V)
1586 J = np.zeros((n_points, n_params), dtype=float)
1587
1588 def get_log_model(pv):
1589 """モデル電流の符号付きlog圧縮表現を計算するヘルパー関数。"""
1590 curr = unpack_free_params(pv, base_params, fix_set)
1591 I_calc = model(
1592 V, curr,
1593 forwardIV=args.forwardIV,
1594 reverseIV=args.reverseIV,
1595 dV=args.dv_blend,
1596 temperature=temperature,
1597 )
1598 return signed_log10_current(I_calc, EPS_I)
1599
1600 f0 = get_log_model(p_free) # 現在のパラメータでのモデル出力
1601
1602 # 各自由パラメータに対して数値微分を実行
1603 for i in range(n_params):
1604 pp = np.copy(p_free)
1605 pp[i] += eps # i番目のパラメータを微小量だけ変更
1606 fi = get_log_model(pp)
1607 J[:, i] = (fi - f0) / eps # 偏微分を近似
1608
1609 return J
1610
1611
1612def estimate_errors(res, V, I_meas, temperature, args, base_params, fix_set):
1613 """
1614 概要: フィッティング結果からパラメータの誤差とモデル電流の信頼区間を推定します。
1615
1616 詳細説明:
1617 ヤコビ行列と残差平方和を用いて共分散行列を計算し、誤差伝播によって
1618 各パラメータの誤差とモデル電流の log10 スケールでの標準偏差を算出します。
1619
1620 引数:
1621 :param res: scipy.optimize.Minimize の結果オブジェクト。最適化されたパラメータ (res.x) と目的関数の値 (res.fun) を含む。
1622 :type res: scipy.optimize.OptimizeResult
1623 :param V: 測定された電圧データ点の配列。
1624 :type V: numpy.ndarray
1625 :param I_meas: 測定された電流データ点の配列。
1626 :type I_meas: numpy.ndarray
1627 :param temperature: シミュレーションの温度 (K)。
1628 :type temperature: float
1629 :param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト。
1630 :type args: argparse.Namespace
1631 :param base_params: 固定パラメータの基準値を含む辞書。
1632 :type base_params: dict
1633 :param fix_set: 固定されたパラメータ名のセット。
1634 :type fix_set: set
1635
1636 戻り値:
1637 :returns: errors_dict (各パラメータ名とその推定誤差の辞書), sigma_log_I (モデル電流のlog10値の標準偏差の配列)
1638 :rtype: tuple[dict, numpy.ndarray]
1639 """
1640 p_opt = np.asarray(res.x, dtype=float) # 最適化された自由パラメータ
1641 n = len(I_meas) # データ点数
1642 p = len(p_opt) # 自由パラメータ数
1643
1644 if p == 0 or n <= p: # 自由パラメータがない、またはデータ点数が少なすぎる場合
1645 return {}, np.zeros(len(V), dtype=float)
1646
1647 J = get_jacobian(p_opt, V, temperature, args, base_params, fix_set) # ヤコビ行列
1648 sigma2_hat = float(res.fun) / max(n - p, 1) # 残差平方和から残差分散を推定
1649
1650 try:
1651 # 共分散行列の計算: (J.T @ J)^-1 * sigma_hat^2
1652 # np.linalg.pinv は特異行列の場合も対応できる擬似逆行列
1653 cov = sigma2_hat * np.linalg.pinv(J.T @ J)
1654 # 自由パラメータ自体の誤差（logスケールの場合）
1655 p_errors_internal = np.sqrt(np.maximum(np.diag(cov), 0.0))
1656 # モデル電流のlog10値の標準偏差
1657 sigma_log_I = np.sqrt(np.maximum(np.diag(J @ cov @ J.T), 0.0))
1658
1659 errors_dict = {}
1660 idx = 0
1661 for name in PARAM_NAMES:
1662 if name in fix_set:
1663 continue
1664 if name in LOG_PARAMS:
1665 # log10変換されたパラメータの誤差を元のスケールに戻す
1666 val = 10.0 ** p_opt[idx]
1667 errors_dict[name] = val * np.log(10.0) * p_errors_internal[idx]
1668 else:
1669 errors_dict[name] = p_errors_internal[idx]
1670 idx += 1
1671
1672 return errors_dict, sigma_log_I
1673
1674 except Exception as e:
1675 print(f"Warning: Error estimation failed ({e})")
1676 return {}, np.zeros(len(V), dtype=float)
1677
1678
1679# =============================================================================
1680# plotting
1681# =============================================================================
1682
1683def initialize_plot(V, I_meas, I_sim, title = None):
1684 """
1685 概要: 最適化プロセスのための初期プロットを設定し、アニメーションオブジェクトを返します。
1686
1687 詳細説明:
1688 リアルタイムでのフィッティング進捗状況を表示するために使用されます。
1689 線形スケールと対数絶対値スケールの2つのサブプロットが作成されます。
1690
1691 引数:
1692 :param V: 電圧データの配列。
1693 :type V: numpy.ndarray
1694 :param I_meas: 測定された電流データの配列。
1695 :type I_meas: numpy.ndarray
1696 :param I_sim: 初期シミュレーション電流データの配列。
1697 :type I_sim: numpy.ndarray
1698 :param title: プロットのメインタイトル。Noneの場合、タイトルは設定されない。
1699 :type title: str or None
1700
1701 戻り値:
1702 :returns: fig (matplotlib.figure.Figure), axes (matplotlib.axes.Axesの配列), anim (アニメーションに使うオブジェクトの辞書)
1703 :rtype: tuple[matplotlib.figure.Figure, numpy.ndarray, dict]
1704 """
1705 anim = {}
1706
1707 plt.ion() # インタラクティブモードをオン
1708 fig, axes = plt.subplots(1, 2, figsize=(11, 5))
1709 ax1, ax2 = axes
1710 if title: ax1.set_title(title)
1711
1712 ax1.plot(V, I_meas, "o", color="black", alpha=0.3, markersize=4, label="data")
1713 line0, = ax1.plot(V, I_sim, "-", color="blue", lw=0.5, label="initial")
1714 line1, = ax1.plot(V, I_sim, "-", color="red", lw=1, label="sim")
1715 band1 = ax1.fill_between(V, I_sim, I_sim, color="lightblue", alpha=0.5) # 初期状態では幅ゼロ
1716 ax1.axhline(0.0, color="k", lw=0.5)
1717 ax1.axvline(0.0, color="k", lw=0.5)
1718 ax1.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
1719 ax1.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
1720 ax1.set_xlabel("V / V", fontsize = fontsize)
1721 ax1.set_ylabel("I / A", fontsize = fontsize)
1722 ax1.grid(True)
1723 ax1.legend()
1724
1725 ax2.set_title("Log abs(I)")
1726 ax2.semilogy(V, np.abs(I_meas) + EPS_I, "o", color="black", alpha=0.3, markersize=4)
1727 line02, = ax2.semilogy(V, np.abs(I_sim) + EPS_I, "-", color="blue", lw=0.5)
1728 line2, = ax2.semilogy(V, np.abs(I_sim) + EPS_I, "-", color="red", lw=1)
1729 band2 = ax2.fill_between(V, np.abs(I_sim) + EPS_I, np.abs(I_sim) + EPS_I, color="lightblue", alpha=0.5) # 初期状態では幅ゼロ
1730
1731 ax2.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
1732 ax2.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
1733 ax2.set_xlabel("V / V", fontsize = fontsize)
1734 ax2.set_ylabel("I / A", fontsize = fontsize)
1735 ax2.grid(True, which="both")
1736
1737 plt.tight_layout()
1738 anim = {
1739 "fig": fig,
1740 "ax1": ax1,
1741 "ax2": ax2,
1742 "line1": line1,
1743 "line2": line2,
1744 "band1": band1,
1745 "band2": band2,
1746 }
1747 return fig, axes, anim
1748
1749def plot_iv(fig, axes, V, I_meas, I_start, I_final, info, mode, sigma_log=None, mode_label = "", outfile=None, pause=False):
1750 """
1751 概要: 測定データとモデルIVカーブをプロットします。
1752
1753 詳細説明:
1754 線形スケールと対数絶対値スケールの2つのサブプロットを作成します。
1755 sigma_log が提供された場合、モデル電流の信頼区間を水色領域で表示します。
1756
1757 引数:
1758 :param fig: プロットに使うmatplotlib.figure.Figureオブジェクト。Noneの場合、新しいFigureを作成。
1759 :type fig: matplotlib.figure.Figure or None
1760 :param axes: プロットに使うmatplotlib.axes.Axesオブジェクトの配列。Noneの場合、新しいAxesを作成。
1761 :type axes: numpy.ndarray or None
1762 :param V: プロットする電圧データの配列。
1763 :type V: numpy.ndarray
1764 :param I_meas: 測定された電流データの配列。Noneの場合、測定データはプロットされない。
1765 :type I_meas: numpy.ndarray or None
1766 :param I_start: 初期モデルによって計算された電流データの配列。Noneの場合、初期モデルはプロットされない。
1767 :type I_start: numpy.ndarray or None
1768 :param I_final: 最終モデルによって計算された電流データの配列。
1769 :type I_final: numpy.ndarray
1770 :param info: ファイル名やレコード時間などの情報を含む辞書。プロットのタイトルに使用。
1771 :type info: dict
1772 :param mode: 現在の実行モード（例: init, fit, sim）を示す文字列。プロットのタイトルに使用。
1773 :type mode: str
1774 :param sigma_log: モデル電流の log10 値の標準偏差の配列。信頼区間を表示するために使用。Noneの場合、信頼区間は表示されない。
1775 :type sigma_log: numpy.ndarray or None
1776 :param mode_label: プロットのタイトルに追加するモードのラベル。
1777 :type mode_label: str
1778 :param outfile: プロットを保存するファイルパス。Noneの場合、保存しない。
1779 :type outfile: str or None
1780 :param pause: プロットを表示後に一時停止するかどうか。Trueの場合、ユーザーのキー入力を待つ。
1781 :type pause: bool
1782
1783 戻り値:
1784 :returns: なし
1785 :rtype: None
1786 """
1787
1788 if fig is None:
1789 fig, axes = plt.subplots(1, 2, figsize=(11, 5))
1790
1791 ax = axes[0]
1792 ax.set_title(f"{info['FileName']} ({mode_label})")
1793 if I_meas is not None:
1794 ax.plot(V, I_meas, "o", color="black", alpha=0.3, markersize=4, label="data")
1795 if I_start is not None:
1796 ax.plot(V, I_start, "-", color="blue", lw=0.5, label="initial")
1797 ax.plot(V, I_final, "-", color="red", lw=1, label="final")
1798
1799 if sigma_log is not None:
1800 log_abs = np.log10(np.abs(I_final) + EPS_I)
1801 upper = 10.0 ** (log_abs + sigma_log)
1802 lower = 10.0 ** (log_abs - sigma_log)
1803 sign = np.sign(I_final) # 電流の符号を保持
1804 ax.fill_between(V, sign * lower, sign * upper, color="lightblue", alpha=0.5)
1805
1806 ax.axhline(0.0, color="k", lw=0.5)
1807 ax.axvline(0.0, color="k", lw=0.5)
1808 ax.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
1809 ax.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
1810 ax.set_xlabel("V / V", fontsize = fontsize)
1811 ax.set_ylabel("I / A", fontsize = fontsize)
1812 ax.grid(True)
1813 ax.legend()
1814
1815 ax = axes[1]
1816 ax.set_title("Log abs(I)")
1817 if I_meas is not None:
1818 ax.semilogy(V, np.abs(I_meas) + EPS_I, "o", color="black", alpha=0.3, markersize=4)
1819 ax.semilogy(V, np.abs(I_final) + EPS_I, "-", color="red", lw=1)
1820
1821 if sigma_log is not None:
1822 log_abs = np.log10(np.abs(I_final) + EPS_I)
1823 upper = 10.0 ** (log_abs + sigma_log)
1824 lower = 10.0 ** (log_abs - sigma_log)
1825 ax.fill_between(V, lower, upper, color="lightblue", alpha=0.5)
1826
1827 ax.tick_params(axis='x', labelsize=fontsize) # x軸の目盛りフォントサイズ
1828 ax.tick_params(axis='y', labelsize=fontsize) # y軸の目盛りフォントサイズ
1829 ax.set_xlabel("V / V", fontsize = fontsize)
1830 ax.set_ylabel("I / A", fontsize = fontsize)
1831 ax.grid(True, which="both")
1832
1833 fig.tight_layout()
1834 if outfile:
1835 print(f"[I/O] Save plot: {outfile}")
1836 fig.savefig(outfile, dpi=160)
1837 if pause:
1838 plt.show()
1839 else:
1840 plt.pause(0.01)
1841
1842
1843def setup_parameters(params, args, csv_fix, args_fix, forwardIV, reverseIV):
1844 """
1845 概要: パラメータ辞書を最終的に設定し、固定されるパラメータを決定します。
1846
1847 詳細説明:
1848 コマンドライン引数で上書きされた値、CSVファイルから読み込まれた値、
1849 デフォルト値、および未使用のモデルメカニズムに基づいて自動的に固定されるパラメータを統合します。
1850
1851 引数:
1852 :param params: 現在のパラメータ値を含む辞書。
1853 :type params: dict
1854 :param args: 解析されたコマンドライン引数を含むオブジェクト。
1855 :type args: argparse.Namespace
1856 :param csv_fix: CSVファイルで固定されたパラメータ名のセット。
1857 :type csv_fix: set
1858 :param args_fix: コマンドライン引数で固定されたパラメータ名のリスト。
1859 :type args_fix: list
1860 :param forwardIV: フォワードバイアスモデルのメカニズムを示す文字列。
1861 :type forwardIV: str
1862 :param reverseIV: リバースバイアスモデルのメカニズムを示す文字列。
1863 :type reverseIV: str
1864
1865 戻り値:
1866 :returns: params (最終的なパラメータ値の辞書), fix_set (最終的に固定されるパラメータ名のセット), used (使用されるパラメータ名のセット), auto_fixed (自動的に固定されたパラメータ名のセット)
1867 :rtype: tuple[dict, set, set, set]
1868
1869 例外:
1870 :raises ValueError: 不明な固定パラメータ名が指定された場合。
1871 """
1872 # コマンドライン引数によるパラメータ上書き
1873 for name in PARAM_NAMES:
1874 val = getattr(args, name)
1875 if val is not None:
1876 params[name] = val
1877
1878 # 不足しているパラメータにデフォルト値を設定
1879 defaults = {
1880 "I0": 1.0e-10,
1881 "ndiode": 1.5,
1882 "IPV": 1.0e-6,
1883 "Rs": 1.0,
1884 "Rsh": 1.0e9,
1885 "A_tfe": 1.0e-30,
1886 "E00": 1.0e-3,
1887 "A_fn": 1.0e-30,
1888 "B_fn": 10.0,
1889 "K_sclc": 1.0e-30,
1890 }
1891 for k, v in defaults.items():
1892 params.setdefault(k, v)
1893
1894 # 固定パラメータの決定
1895 fix_set = set(csv_fix) | set(args_fix) # CSVとコマンドライン引数の両方で指定されたものを結合
1896
1897 # 使用されていないモデルパラメータを自動的に固定
1898 used = used_params_from_modes(forwardIV, reverseIV)
1899 auto_fixed = set(PARAM_NAMES) - used
1900 fix_set |= auto_fixed
1901
1902 # 不明な固定パラメータ名のチェック
1903 unknown_fix = [x for x in fix_set if x not in PARAM_NAMES]
1904 if unknown_fix:
1905 raise ValueError(f"Unknown parameter name(s) in --fix: {unknown_fix}")
1906
1907 print(f"Used parameters : {' '.join(sorted(used))}")
1908 print(f"Auto-fixed unused : {' '.join(sorted(auto_fixed)) if auto_fixed else '(none)'}")
1909
1910 return params, fix_set, used, auto_fixed
1911
1912
1913def print_parameter_repeat(params, fix_set, used, auto_fixed, title=""):
1914 """
1915 概要: 現在のモデルパラメータとその固定状態をコンソールに出力します。
1916
1917 詳細説明:
1918 これは、print_args_and_derived の「Derived values」セクションに似ていますが、
1919 最終的に決定された各パラメータの具体的な値と、それが固定されているかどうかの状態を詳細に示します。
1920
1921 引数:
1922 :param params: 最終的なパラメータ値を含む辞書。
1923 :type params: dict
1924 :param fix_set: 固定されるパラメータ名のセット。
1925 :type fix_set: set
1926 :param used: モデルで使用されるパラメータ名のセット。
1927 :type used: set
1928 :param auto_fixed: モデルによって自動的に固定されたパラメータ名のセット。
1929 :type auto_fixed: set
1930 :param title: 出力の前に表示されるタイトル文字列。
1931 :type title: str
1932
1933 戻り値:
1934 :returns: なし
1935 :rtype: None
1936 """
1937 if title:
1938 print(title)
1939 for name in PARAM_NAMES:
1940 status = ""
1941 if name in fix_set:
1942 status = "(Fixed)"
1943 if name in auto_fixed:
1944 status = "(Auto-Fixed)"
1945
1946 print(f" {name:8s}: {params.get(name, 'N/A'):.6e} {status}")
1947 print()
1948
1949
1950# =============================================================================
1951# fit execution
1952# =============================================================================
1953
1954def exec_fit(V, I_meas, info, csv_path, params, fix_set, args):
1955 """
1956 概要: fit モードの本体処理。最適化、誤差推定、保存、最終プロットまでをまとめて実行する。
1957
1958 詳細説明:
1959 この関数は、scipy.optimize.minimize を用いてSDMパラメータを測定データにフィッティングします。
1960 フィッティング中には、コンソールへの進捗表示とリアルタイムのIVカーブアニメーションが行われます。
1961 最適化後、パラメータの誤差が推定され、結果がCSVファイルに保存され、最終的なプロットが表示されます。
1962
1963 引数:
1964 :param V: 測定された電圧データ点の配列。
1965 :type V: numpy.ndarray
1966 :param I_meas: 測定された電流データ点の配列。
1967 :type I_meas: numpy.ndarray
1968 :param info: ファイル名やレコード時間などの情報を含む辞書。プロットのタイトルに使用。
1969 :type info: dict
1970 :param csv_path: パラメータを保存するCSVファイルのパス。
1971 :type csv_path: pathlib.Path
1972 :param params: 初期パラメータ値を含む辞書。
1973 :type params: dict
1974 :param fix_set: 固定されるパラメータ名のセット。
1975 :type fix_set: set
1976 :param args: コマンドライン引数パーサーから得られた設定を含むオブジェクト。
1977 :type args: argparse.Namespace
1978
1979 戻り値:
1980 :returns: params (最適化後のパラメータ値の辞書), rss (残差平方和), errors (各パラメータの推定誤差の辞書), sigma_log_curve (モデル電流のlog10値の標準偏差の配列)
1981 :rtype: tuple[dict, float, dict, numpy.ndarray]
1982 """
1983
1984 print(f"Execute {args.mode} mode.")
1985 params_init = dict(params)
1986 I_start = model(V, params, args.forwardIV, args.reverseIV, args.dv_blend, args.temperature)
1987 fig, axes, anim = initialize_plot(V, I_meas, I_start, title="Optimization progress")
1988
1989 p0 = pack_free_params(params, fix_set)
1990 state = {"count": 0}
1991
1992 def callback(xk):
1993 state["count"] += 1
1994 curr = unpack_free_params(xk, params, fix_set)
1995
1996 if state["count"] % max(args.ninterval_print, 1) == 0:
1997 rss = objective(xk, V, I_meas, args, params, fix_set)
1998 pstr = " ".join(
1999 f"{n}={curr[n]:.3e}" for n in PARAM_NAMES if n not in fix_set
2000 )
2001 print(f"{state['count']:4d}: RSS={rss:.6e} {pstr}")
2002
2003 if state["count"] % max(args.ninterval_plot, 1) == 0:
2004 I_anim = model(
2005 V, curr,
2006 forwardIV=args.forwardIV,
2007 reverseIV=args.reverseIV,
2008 dV=args.dv_blend,
2009 temperature=args.temperature,
2010 )
2011 anim["line1"].set_ydata(I_anim)
2012 anim["line2"].set_ydata(np.abs(I_anim) + EPS_I)
2013 anim["ax1"].relim()
2014 anim["ax1"].autoscale_view()
2015 anim["ax2"].relim()
2016 anim["ax2"].autoscale_view()
2017 anim["fig"].canvas.draw()
2018 anim["fig"].canvas.flush_events()
2019 plt.pause(0.001)
2020
2021 if len(p0) == 0:
2022 print("All parameters are fixed. Skipping optimization.")
2023 class DummyResult:
2024 x = np.array([], dtype=float)
2025 fun = objective(np.array([], dtype=float), V, I_meas, args, params, fix_set)
2026 res = DummyResult()
2027 rss = float(res.fun)
2028 else:
2029 res = minimize(
2030 objective,
2031 p0,
2032 args=(V, I_meas, args, params, fix_set),
2033 method=args.method,
2034 callback=callback,
2035 options={
2036 "maxiter": int(args.nmaxiter),
2037 "xatol": float(args.tol),
2038 "fatol": float(args.tol),
2039 "disp": False,
2040 },
2041 )
2042 params = unpack_free_params(res.x, params, fix_set)
2043 rss = float(res.fun)
2044
2045 errors, sigma_log_curve = estimate_errors(
2046 res, V, I_meas, args.temperature, args, params, fix_set
2047 )
2048
2049 save_param_csv(csv_path, params, fix_set, errors=errors, rss=rss)
2050
2051 plt.ioff()
2052
2053 print("\n--- Fit Results ---")
2054 print(f"RSS_logy : {rss:.6e}")
2055 for n in PARAM_NAMES:
2056 if n in errors:
2057 print(f" {n:8s}: {params[n]:.6e} +/- {errors[n]:.2e}")
2058 else:
2059 print(f" {n:8s}: {params[n]:.6e}" + (" (Fixed)" if n in fix_set else ""))
2060
2061 print()
2062 print("Parameters estimated from input IV data (Use if the fitting accurary is poor)")
2063 est = estimate_initial_params(V, I_meas, args.temperature)
2064 print(f" Rs_est : {est['Rs']:.6f} ohm")
2065 print(f" Rsh_est : {est['Rsh']:.6f} ohm")
2066 print(f" ISC_est : {est['ISC_est']:.6e} A")
2067 print(f" VOC_est : {est['VOC_est']:.6f} V")
2068
2069 I_final = model(
2070 V, params,
2071 forwardIV=args.forwardIV,
2072 reverseIV=args.reverseIV,
2073 dV=args.dv_blend,
2074 temperature=args.temperature,
2075 )
2076
2077 # Excel保存
2078 outfile_xlsx = f"{args.outprefix}-fitted.xlsx"
2079 save_iv_to_excel(
2080 outfile_xlsx,
2081 V=V,
2082 I_input=I_meas,
2083 I_init=I_start,
2084 I_fit=I_final,
2085 info=info,
2086 params_init=params_init,
2087 params_final=params,
2088 rss=rss,
2089 errors=errors,
2090 )
2091
2092 # 画像保存してから表示
2093 outfile_iv = f"{args.outprefix}-iv.png"
2094 plot_iv(
2095 fig, axes, V, I_meas, I_start, I_final, info, args.mode,
2096 sigma_log=sigma_log_curve,
2097 outfile=outfile_iv,
2098 pause=False
2099 )
2100 input("\nPress ENTER to terminate\n")
2101
2102 return params, rss, errors, sigma_log_curve
2103
2104
2105# =============================================================================
2106# main
2107# =============================================================================
2108
2109_original_print = None
2110_redirect_fp = None
2111def dual_print(*args, **kwargs):
2112 """
2113 概要: 標準出力とログファイルの両方にメッセージを出力します。
2114
2115 詳細説明:
2116 組み込みの print 関数をオーバーライドするために使用され、
2117 すべての出力がコンソールと指定されたログファイルの両方に書き込まれます。
2118
2119 引数:
2120 :param args: print関数に渡される可変長引数。
2121 :type args: tuple
2122 :param kwargs: print関数に渡されるキーワード引数。
2123 :type kwargs: dict
2124
2125 戻り値:
2126 :returns: なし
2127 :rtype: None
2128 """
2129 _original_print(*args, **kwargs)
2130
2131 file_kwargs = dict(kwargs)
2132 file_kwargs.pop("file", None) # file指定があっても無視してログファイルへ
2133 file_kwargs.setdefault("flush", True)
2134
2135 _original_print(*args, file=_redirect_fp, **file_kwargs)
2136
2137
2138def main():
2139 """
2140 概要: スクリプトのエントリーポイント。コマンドライン引数を解析し、IVフィッティングツールを実行します。
2141
2142 詳細説明:
2143 この関数は、初期値推定 (init)、モデルフィッティング (fit)、またはシミュレーション (sim) の
2144 いずれかのモードで動作します。データ読み込み、パラメータ管理、最適化、結果のプロットを行います。
2145
2146 引数:
2147 なし
2148
2149 戻り値:
2150 :returns: なし
2151 :rtype: None
2152 """
2153
2154 global _original_print, _redirect_fp
2155
2156 print()
2157 print("Read commandline arguments")
2158 args, parser = initialize()
2159
2160 validate_args(args)
2161
2162 _original_print = builtins.print
2163 logfile = f"{args.outprefix}_{args.mode}.txt"
2164 print(f"Open log file [{logfile}]")
2165 _redirect_fp = open(logfile, "w", encoding="utf-8")
2166 builtins.print = dual_print
2167
2168 print_args_and_derived(args)
2169 print(f"forwardIV model : {args.forwardIV}")
2170 print(f"reverseIV model : {args.reverseIV}")
2171
2172 # データ入力処理
2173 print()
2174 if args.infile:
2175 print(f"Read data from [{args.infile}]")
2176 print(f" x range: {args.xmin} - {args.xmax}, every {args.ndataskip} data")
2177 V, I_meas, info = read_data(args.infile, xmin=args.xmin, xmax=args.xmax, ndataskip=args.ndataskip)
2178 csv_path = Path(args.infile).stem + "-parameters.csv"
2179 else:
2180 if args.mode != "sim":
2181 raise ValueError("--infile is required for mode init/fit/analyze.")
2182
2183 print(f"Create simulation data")
2184 # シミュレーションモードの場合、ダミーデータを生成
2185 V = np.linspace(-0.2, 0.8, 500)
2186 I_meas = None
2187 info = {"FileName": "Simulation", "RecordTime": "N/A"}
2188 csv_path = "sim-parameters.csv"
2189
2190 # パラメータの読み込みと初期値設定
2191 print()
2192 if os.path.exists(csv_path):
2193 print(f"Read parameters from [{csv_path}]")
2194 params, csv_fix = load_param_csv(csv_path)
2195 else:
2196 print(f"Parameter file [{csv_path}] does not exist. Use default.")
2197 params = {}
2198 csv_fix = set()
2199
2200 if (args.mode == "init" or args.mode == "analyze") or not any(name in params for name in PARAM_NAMES):
2201 if I_meas is not None:
2202 print()
2203 print(f"Estimate initial parameters at T={args.temperature}")
2204 est = estimate_initial_params(V, I_meas, args.temperature)
2205 else:
2206 print()
2207 print(f"Use default initial parameters at T={args.temperature}")
2208 est = DEFAULT_PARAMETERS
2209 params.update(est)
2210
2211 print()
2212 print("Final setup model parameters from paramter file and commandline arguments")
2213 params, fix_set, used, auto_fixed = setup_parameters(params, args, csv_fix, args.fix, args.forwardIV, args.reverseIV)
2214 print_parameter_repeat(params, fix_set, used, auto_fixed, title="=== Final model parameters ===")
2215
2216 print()
2217 if args.mode == "analyze" or args.mode == "init":
2218 print(f"Execute {args.mode} mode.")
2219 print(" Extracted parameters")
2220 print(f" ISC_est : {params.get('ISC_est', np.nan):.6e} A")
2221 print(f" VOC_est : {params.get('VOC_est', np.nan):.6f} V")
2222 print(f" Vsh_est : {params.get('Vsh_est', np.nan):.6f} V")
2223 print(f" Ish_est : {params.get('Ish_est', np.nan):.6e} A")
2224 print(f" Vtfe_est: {params.get('Vtfe_est', np.nan):.6f} V")
2225 print(f" Itfe_est: {params.get('Itfe_est', np.nan):.6e} A")
2226
2227 I_model = model(
2228 V, params,
2229 forwardIV=args.forwardIV,
2230 reverseIV=args.reverseIV,
2231 dV=args.dv_blend,
2232 temperature=args.temperature,
2233 )
2234
2235 if args.mode == "init":
2236 save_param_csv(csv_path, params, fix_set)
2237 outfile_xlsx = f"{args.outprefix}-fitted.xlsx"
2238 save_iv_to_excel(
2239 outfile_xlsx,
2240 V=V,
2241 I_input=I_meas,
2242 I_init=None,
2243 I_fit=I_model,
2244 info=info,
2245 params_init=params,
2246 params_final=None,
2247 rss=None,
2248 errors=None,
2249 )
2250
2251 plot_iv(
2252 None, None, V, I_meas, None, I_model, info, args.mode,
2253 sigma_log=None,
2254 outfile=f"{args.outprefix}-iv.png",
2255 pause=False,
2256 )
2257 input("\nPress ENTER to terminate\n")
2258
2259 elif args.mode == "fit":
2260 params, rss, errors, sigma_log_curve = exec_fit(
2261 V, I_meas, info, csv_path, params, fix_set, args
2262 )
2263
2264 elif args.mode == "sim":
2265 print(f"Execute {args.mode} mode.")
2266 I_model = model(
2267 V, params,
2268 forwardIV=args.forwardIV,
2269 reverseIV=args.reverseIV,
2270 dV=args.dv_blend,
2271 temperature=args.temperature,
2272 )
2273
2274 outfile_xlsx = f"{args.outprefix}-fitted.xlsx"
2275 save_iv_to_excel(
2276 outfile_xlsx,
2277 V=V,
2278 I_input=None,
2279 I_init=None,
2280 I_fit=I_model,
2281 info=info,
2282 params_init=params,
2283 params_final=None,
2284 rss=None,
2285 errors=None,
2286 )
2287
2288 plot_iv(
2289 None, None, V, I_meas, None, I_model, info, args.mode,
2290 sigma_log=None,
2291 outfile=f"{args.outprefix}-iv.png",
2292 pause=False,
2293 )
2294 input("\nPress ENTER to terminate\n")
2295
2296
2297if __name__ == "__main__":
2298 try:
2299 main()
2300 except Exception:
2301 traceback.print_exc()
2302 sys.exit(1)
2303 finally:
2304 if _redirect_fp: _redirect_fp.close()