next up previous index
Next: B.5 SPICEのモデル Up: B. 真空管のモデル Previous: B.3 Rによるインプリメント


B.4 パラメータの求め方(キャリブレーション)

この節では,プレート特性図から読み取ったデータから, 前節のモデルのパラメータを推定する方法を説明します.

B.4.1 データのフォーマット

プレート特性図から何点かのデータを読み取り, ファイルに格納しておきます. ファイルの形式はCSVで, グリッド電圧,プレート電圧,プレート電流の3列のデータです. 1行目は見出しでデータとしては扱われませんが, 例のように Eg,Ep,Ip としておいてください. 単位は,VおよびAで入力します.

ファイルの例を次に示します.

    1 Eg,Ep,Ip
    2 0,0,0
    3 0,50,0.005
    4 0,100,0.0118
    5 0,150,0.02
    6 -2,17,0
    7 -2,50,0.0015
    8 -2,100,0.0062
    9 -2,150,0.013
   10 -2,200,0.021
   11 -4,43,0
   12 -4,100,0.0027
   13 -4,150,0.0075
   14 -4,200,0.0147
   15 -4,250,0.0232
   16 -6,70,0
   17 -6,100,8.0E-4
   18 -6,150,0.0040
   19 -6,200,0.0094
   20 -6,250,0.017
   21 -8,95,0
   22 -8,150,0.0017
   23 -8,200,0.0053
   24 -8,250,0.0117
   25 -8,300,0.02
   26 -10,125,0
   27 -10,150,4.0E-4
   28 -10,200,0.0027
   29 -10,250,0.0074
   30 -10,300,0.0142
   31 -12,157,0
   32 -12,200,0.0012
   33 -12,250,0.0043
   34 -12,300,0.0097
   35 -12,350,0.0164
   36 -14,183,0
   37 -14,200,3.0E-4
   38 -14,250,0.0024
   39 -14,300,0.0061
   40 -14,350,0.0115
   41 -16,205,0
   42 -16,250,0.0010
   43 -16,300,0.0034
   44 -16,350,0.0077
   45 -18,228,0
   46 -18,250,3.0E-4
   47 -18,300,0.0017
   48 -18,350,0.0048
   49 -20,265,0
   50 -20,300,7.0E-4
   51 -20,350,0.0028
この例では,50点のプレート電流を読み取っていますが, 20点程度でも十分な精度のパラメータを得ることができます. 各グリッド電圧に対して,最低でもカットオフに近いプレート電圧と, それよりもかなり大きな電圧(特性図上に載っている最大のプレート電圧あるいは プレート損失近辺のプレート電圧)におけるプレート電流を入力しておきます.

CSVファイルは,テキストエディタで作成することができます. または,Microsoft Excelなどの表計算ソフトで作成し, CSV形式で保存することもできます. 表計算ソフトを用いれば,プレート電流をmA単位で入力して, A単位に変換することも容易でしょう.

B.4.2 キャリブレーションを行なう関数 Ip.cal

    1 "Ip.cal" <-
    2 function(fn, name, gfn, G.lim=NULL, Ig.ratio=NULL, verbose=FALSE, ...)
    3 {
    4     # モデルのパラメータを求める
    5     # fn: データファイル名
    6     # name: パラメータオブジェクト名
    7     # gfn: キャリブレーションの状況を示すプレート特性のグラフファイル名.
    8     #   指定されていなければグラフを画面に表示する.
    9     # G.lim: Ep=Eg のパービアンス(省略可)
   10     # Ig.ratio: Ep=Eg のときのカソード電流に対するグリッド電流の比率(省略可)
   11     # verbose: 途中経過の表示
   12     # ...: 電極間容量および他のパラメータ
   13     #   Cgp, Cgk Cpk, Ci, Co: 電極間容量
   14     #   ig=(Ep,Ip,Ig): グリッド電流データ
   15     #   pentode.chara=(Ep,Eg2,Eg,Ip,Ig2,rp): 五極管動作例
   16 
   17     sig <- function(x) signif(x, 8)
   18     print.level <- 1        # 診断出力のレベル
   19 
   20     # キャリブレーション
   21     x <- scan(fn, sep=",", what="")     # データファイルを読む
   22     egepip <- matrix(as.numeric(x[-(1:3)]), ncol=3, byrow=TRUE)
   23     eg <- egepip[,1]        # グリッド電圧
   24     ep <- egepip[,2]        # プレート電圧
   25     ip <- egepip[,3]        # プレート電流
   26     n <- length(eg)         # キャリブレーションデータのポイント数
   27 
   28     if (!missing(gfn))
   29         psi(file=gfn, wh=c(103.9, 68.8))
   30     
   31     # グラフを描く範囲を決定する
   32     xlim <- c(0, max(ep)*1.04)
   33     ylim <- c(0, max(ip)*1.04)
   34     xlab <- eunit1(pretty(xlim, 8), "V")
   35     ylab <- eunit1(pretty(ylim), "A")
   36 
   37     # キャリブレーションポイントをグラフに描く(ここではまだ描かない)
   38     plot(ep, ip, type="n", xlim=xlim, ylim=ylim,
   39         xaxs="i", yaxs="i", xaxt="n", yaxt="n",
   40         xlab=paste("Ep", xlab$unit),
   41         ylab=paste("Ip", ylab$unit))
   42 
   43     # 同じグリッド電圧のポイントを同じ色で描く
   44     ueg <- unique(eg)               # 一意のグリッド電圧
   45     for (i in 1:length(ueg)) {
   46         ok <- eg == ueg[i]
   47         points(ep[ok], ip[ok], pch="+", col=(i-1)%%4+3)
   48     }
   49 
   50     # 軸を描く
   51     axis(1, at=xlab$at, lab=xlab$lab)
   52     axis(1, at=xlab$at, lab=FALSE, tck=1, lty=2, lwd=0.5)
   53     axis(2, at=ylab$at, lab=ylab$lab)
   54     axis(2, at=ylab$at, lab=FALSE, tck=1, lty=2, lwd=0.5)
   55 
   56     if (is.null(G.lim)) {
   57         # 最適化の変数からモデルの変数への変換関数
   58         decode <- function(x) {
   59             z <- exp(x)
   60             z[3] <- z[3] / (z[3] + 1) * 0.65 + 0.35
   61             z[4] <- z[4] / (z[4] + 1) * 2 - 1
   62             list(G=z[1], muc=z[2], alpha=z[3], Ego=z[4])
   63         }
   64     } else {
   65         # 最適化の変数からモデルの変数への変換関数
   66         decode <- function(x) {
   67             z <- exp(x)
   68             z[3] <- z[3] / (z[3] + 1) * 0.65 + 0.35
   69             z[4] <- z[4] / (z[4] + 1) * 2 - 1
   70             list(G=z[1], muc=z[2], alpha=z[3], Ego=z[4], G.lim=G.lim, Ig.ratio=Ig.ratio)
   71         }
   72     }
   73 
   74     # モデルの変数から最適化の変数への変換関数
   75     encode <- function(x) {
   76         x[3] <- (x[3] - 0.35) / 0.65
   77         x[3] <- x[3] / (1 - x[3])
   78         x[4] <- (x[4] + 1) / 2
   79         x[4] <- x[4] / (1 - x[4])
   80         log(x)
   81     }
   82 
   83     # 最適化する関数
   84     f <- function(x) {
   85         p <- decode(x)                  # モデルの変数へ変換
   86         if (print.level == 2)
   87             str(p)
   88         ipe <- Ip(p, Eg=eg, Ep=ep)      # モデルによるプレート電流
   89         if (any(is.na(ipe))) {          # モデルからNAが出てきた時の処置
   90             print(egepip[, 1:2])        # (ここにくるはずはない?)
   91             print(ipe)
   92             str(p)
   93             stop()
   94         }
   95         ipd <- ipe^(2/3) - ip^(2/3)         # モデルと実際の値の差
   96         sqrt(mean(ipd^2))/mean(ip^(2/3))    # プレート電流の平均値に対する
   97                                             # RMS偏差の大きさ
   98     }
   99 
  100     # パラメータ初期値の推定
  101     # 増幅率の推定
  102     mineg <- min(eg)                    # データで最も深いグリッド電圧
  103     if (mineg < 0) {
  104         minep <- min(ep[eg == mineg])   # その中で最もプレート電圧が低いもの
  105         muhat <- minep / -mineg         # それをカットオフとみなして mu_c を推定
  106     } else
  107         muhat <- 30                     # B級、ポジティブグリッド管用
  108 
  109     # パービアンスの推定
  110     mu0 <- 5/3 * muhat                  # Eg=0の mu
  111     ep0 <- ifelse(eg == max(eg), ep, 0)
  112     idx <- ep0 == max(ep0)              # Eg=0で最もプレート電圧が高いものにより
  113     Ghat <- 2/3 * ip[idx] / (0.6 + ep[idx]/mu0)^1.5     # パービアンスを推定
  114     if (verbose) {
  115         cat("muhat=", muhat, "\n", sep=" ")
  116         cat("Ghat=", Ghat, "\n", sep=" ")
  117     }
  118 
  119     ini <- encode(c(Ghat, muhat, 0.6, 0.6))     # 他のパラメータの初期値は0.6
  120 
  121     # 最適化
  122     optim.out <- optim(ini, f,
  123                     control=list(parscale=ini, maxit=5000, reltol=1e-10))
  124     cat("code=", optim.out$convergence, " err=", optim.out$value, "\n", sep="")
  125 
  126     # パラメータオブジェクトの作成と保存
  127     est <- decode(optim.out$par)            # 推定されたモデルのパラメータ
  128     for (i in 1:length(est))                # 有効数字を8桁にする
  129         est[[i]] <- sig(est[[i]])
  130     est <- c(est, list(...))
  131     est$err <- optim.out$value
  132     est$code <- optim.out$convergence
  133     if (!is.null(est$pentode.chara)) {
  134         pen.Ep <- est$pentode.chara[1]
  135         pen.Eg2 <- est$pentode.chara[2]
  136         pen.Eg <- est$pentode.chara[3]
  137         pen.Ip <- est$pentode.chara[4]
  138         pen.Ig2 <- est$pentode.chara[5]
  139         pen.rp <- est$pentode.chara[6]
  140         g2r <- pen.Ig2/(pen.Ip + pen.Ig2)
  141         a <- (1-pen.Ep/(pen.Ep+10))^1.5
  142         r <- (g2r - a)/(1 - a)              # (1.41)
  143         est$g2.r <- sig(r)
  144         est$Ea <- sig(pen.Eg2 - pen.Ip * pen.rp * 1.5)
  145     }
  146     if (!missing(name))                     # オブジェクト名が指定されていれば
  147         assign(name, est, env=.GlobalEnv)   # 大域オブジェクトとする
  148 
  149     # モデルによるプレート特性をグラフに書き加える
  150     g.plate.cal(est, ueg, xlim, ylim)
  151 
  152     if (!missing(gfn))
  153         dev.off()
  154     invisible()
  155 }

まず fn で指定されたファイルからデータを読み込み, 各列ごとに分解して eg, ep, ip に付値しておきます.

グラフを描く範囲は,最大のプレート電圧やプレート電流の4%増しとして, すべての点がグラフの枠内に収まるようにしています.

最適化を行なう関数 optim は, パラメータの定義域として - $ \infty$ から $ \infty$ の範囲を探索しますが, モデルのパラメータは0以上なので, 最適化のパラメータとモデルのパラメータの間で 相互に変換する関数を定義しておきます. それが関数 decodeencode で, 最適化のパラメータの指数をとったものがモデルのパラメータ(> 0)となり, モデルのパラメータの対数をとったものが最適化のパラメータとなります. また,モデルのパラメータの対数をとることにより, パラメータ間の桁数の違いもかなり吸収され, パラメータのスケールの調整に気を遣わなくてすみます. パラメータ alpha の範囲は0.35から1であり, パラメータ Ego の範囲は -1 から1なので,変換が異なります.

最適化の対象となる関数は f で, 最適化関数から渡されたベクトルをモデルのパラメータに変換し, モデルによりプレート電流の1.5乗根を求め,実際のプレート電流の1.5乗根との差をとります. この差はプラスにもマイナスにもなりうるので, 誤差の2乗の平均値の平方根をとり, プレート電流の1.5乗根の平均に対する比をできるだけ小さくなるようにします. このように標準化したのは, 異なる真空管に対するフィッティングの良さを比較できるようにするためです. プレート電流の1.5乗根をとるのは, 各グリッド電圧に対するプレート電圧とプレート電流の組ができるだけ直線上に並ぶようにするためです. 1.5乗根をとることにより,カットオフ近辺の誤差が小さくなるように最適化されます.

最適化に入る前に,パラメータの初期値を推定します. $ \alpha$Ego は先験的に値がわかっているので,それぞれ0.6とします. $ \mu_{c}^{}$ については, 最もグリッド電圧が低い場合のカットオフとなるプレート電圧から求めます. またパービアンスについては,Eg = 0 の特性曲線から求めます.

推定されたパラメータは有効数字を8桁に制限しています. これは,この関数を用いてパラメータの推定をせずに, 本書に載っているパラメータの数値を使用しても 同じ結果が再現されるようにするためです.

Ip.cal の使用例を以下に示します. ここでは,12AU7.csv というファイルからデータを読み込み, t12AU7 というオブジェクトにパラメータを作成します. グラフは画面に表示されます.

> Ip.cal("12AU7.csv", "t12AU7", Cgp=1.5e-12, Ci=1.6e-12, Co=0.4e-12)
Read 153 items          # 読み込んだデータの総数(1行目を含む)
code=0 err=0.02995393   # 収束した
> str(t12AU7)
List of 9
 $ G    : num 0.000548
 $ muc  : num 13.2
 $ alpha: num 0.58
 $ Ego  : num 0.829
 $ Cgp  : num 1.5e-12
 $ Ci   : num 1.6e-12
 $ Co   : num 4e-13
 $ err  : num 0.0300
 $ code : int 0

電極間容量は,次のパラメータ名で与えます.

パラメータ名 意味
Cgp グリッド-プレート間容量
Cgk グリッド-カソード間容量
Cpk プレート-カソード間容量
Ci 入力容量
Co 入力容量
パラメータが不足して計算できない場合は 1pF が使用されます.

グリッド電流のデータがある場合は, プレート電圧(=グリッド電圧),プレート電流,グリッド電流の3つの値をベクトルとして パラメータ名 ig で与えます. たとえば 801 の場合, Ep = Eg = 100 V のとき Ip = 225 mA, Ig = 48 mA ですので,

> Ip.cal("801.csv", "t801", Cgp=6e-12, Cgk=4.5e-12, Cpk=1.5e-12,
    ig=c(100, 225e-3, 48e-3))
のようにして与えます.

五極管の場合は,標準的な動作例からプレート電圧,スクリーングリッド電圧,コントロールグリッド電圧,プレート電流,スクリーングリッド電流,プレート内部抵抗の6つの値をベクトルとしてパラメータ名 pentode.chara で与えます. たとえば 6L6 の場合, Ep = 250 V, Eg2 = 250 V, Eg = - 14 V, Ip = 72 mA, Ig2 = 5 mA, rp = 22.5 kΩ ですから,

> Ip.cal("6L6.csv", "t6L6T", Cgp=0.6e-12, Ci=10e-12, Co=6.5e-12,
    pentode.chara=c(250, 250, -14, 72e-3, 5e-3, 22.5e3),
    ig=c(10, 18e-3, 32e-3))
のようにして与えます. できるだけ Ep = Eg2 のデータを与えます. ただし,プレート内部抵抗の値については, プレート特性図ができるだけ再現されるように適宜調整します.


next up previous index
Next: B.5 SPICEのモデル Up: B. 真空管のモデル Previous: B.3 Rによるインプリメント
Ayumi Nakabayashi
平成19年6月28日