Tensorflowの環境を構築して線形回帰を試す

年齢と点数の関係

GoogleやMicrosoft、IBMが提供する機械学習済みのデータを使うのであれば、自身で環境を構築する必要はありません。

しかし、独自の学習をさせるのであれば自身のコンピュータに環境を構築するか、或いはGoogleなどが提供するクラウド上の機械学習の環境が必要です。

以下のサイトを参考にさせていただき、MacOSにTensorflowをインストールして環境を構築しました。
https://qiita.com/junichiro/items/8886f3976fc20f73335f

上記のサイトでは、予め線形関数の係数を定めて、この関数からデータを作成して学習データとしています。

今回は、下記のサイトのデータを使わせていただき、線形回帰の機械学習を試してみました。
https://kyoto-edu.sakura.ne.jp/?&course=statistics&content=correl

実行結果は以下のとおり。

python tensorflow-test.py

0 [8.540184] [0.9793252]
200 [6.1199374] [-0.58089983]
400 [6.234987] [-1.5768614]
600 [6.325655] [-2.3617592]
800 [6.3971086] [-2.9803212]
1000 [6.4534197] [-3.467797]
1200 [6.4977975] [-3.8519664]
1400 [6.5327706] [-4.1547227]
1600 [6.560332] [-4.3933163]
1800 [6.582052] [-4.581347]
2000 [6.5991697] [-4.7295313]

0歳の得点がマイナスというのも変なモデルですが、そこはそういうデータだったということで良しとします。

Pythonのソースファイルは以下です。

import tensorflow as tf
import numpy as np

x_data = np.loadtxt('exam-result.csv', delimiter=',', usecols=[3], skiprows=1)
y_data = np.loadtxt('exam-result.csv', delimiter=',', usecols=[1], skiprows=1)

# Try to find values for W and b that compute y_data = W * x_data + b
# (We know that W should be 0.1 and b 0.3, but TensorFlow will
# figure that out for us.)
W = tf.Variable(tf.zeros([1]))
b = tf.Variable(tf.zeros([1]))
y = W * x_data + b

# Minimize the mean squared errors.
loss = tf.reduce_mean(tf.square(y - y_data))
optimizer = tf.train.GradientDescentOptimizer(0.01)
train = optimizer.minimize(loss)

# Before starting, initialize the variables.  We will 'run' this first.
init = tf.global_variables_initializer()

# Launch the graph.
sess = tf.Session()
sess.run(init)

# Fit the line.
for step in range(2001):
    sess.run(train)
    if step % 200 == 0:
        print(step, sess.run(W), sess.run(b))

# Close the Session when we're done.
sess.close()

AEiスクリプトで埋込SRモデルの解析を連続実行する

動的サブストラクチャ解析の一つに、建屋の一部(または全部)が表層地盤に埋め込まれた構造を、地盤の応答解析の結果を引き継ぎ、建屋の応答解析を実施するという方法があります。

この場合、建屋側面地盤の非線形性を等価線形解析で求め、この物性で求めた建屋側面のインピーダンスから地盤剛性と地盤減衰を算出して、建屋に接続して解析を行うことが一般的です。

以下は、その一連の解析を、入力地震動を変更しながら、また地盤物性と建屋減衰を変更してAEiスクリプトで連続実行するデモです。

入力地震動はエルセントロNSと日本建築センター模擬波(基盤波)、八戸NSの3種、地盤物性は標準、+σ、-σの3種、建屋減衰は5%と3%の2種で合計18ケースの解析を実施します。

なお、地盤の応答解析以外は自前のプログラムを使用しております。

地盤の応答解析は吉田望氏が作成されたDYNEQ336を使用させていただきました。
たいへん使いやすく、ドキュメントもしっかりいて素晴らしいものだと思います。

デモ動画のOSはMacOSですが、AEiスクリプトはWindowsでも、Linuxでも動作します。

参考として以下にaeiスクリプトを記載します。

print " aei スクリプトで、入力地震動と地盤物性を変化させ、埋込SRを連続解析して各ケースの解析結果のレポートを作成する"
pause

print "解析定義のExcelから解析ケースのシートを抽出する"
cmd{
	excelcsv -proc ex2csv -excel analysis.xlsx -sheet case -csv case.csv
	cat case.csv
	}
	pause

# 連続解析のループ
for i = 1 to 18
	L = i + 1

	print "解析ケースのCSVから入力波と地盤の情報を変数に代入する"
	csvfile = case.csv
	$NUM  = csv_r[L]_c[1]
	$SOIL = csv_r[L]_c[2]
	$WAVE = csv_r[L]_c[3]
	$WNUM = csv_r[L]_c[4]
	$WSTP = csv_r[L]_c[5]
	$BLDH = csv_r[L]_c[6]
	print "ケースNo   = ",$NUM
	print "地盤条件   = ",$SOIL
	print "入力地震動 = ",$WAVE,"( ",$WSTP,")"
	print "建屋減衰   = ",$BLDH
	pause

	print "各解析データのファイル名を定義する"
	$CASE = CASE<$NUM>
	$INPUT_WAV = <$WAVE>.wav
	$JBNWAV    = <$SOIL>_<$WAVE>
	$SHAKE_DAT = Shake_<$JBNWAV>.dat
	$SHAKE_OUT = Shake_<$JBNWAV>.out
	$SHAKE_WAV = Shake_<$JBNWAV>.wav
	$NOVAK_DAT = Novak_<$JBNWAV>.dat
	$NOVAK_OUT = Novak_<$JBNWAV>.out
	$ADMIT_DAT = Admit_<$SOIL>.dat
	$ADMIT_OUT = Admit_<$SOIL>.out
	$NAMAZ_DAT = <$CASE>_Namaz_<$JBNWAV>_h<$BLDH>.dat
	$NAMAZ_OUT = <$CASE>_Namaz_<$JBNWAV>_h<$BLDH>.out
	$NAMAZ_WAV = <$CASE>_Namaz_<$JBNWAV>_h<$BLDH>.wav
	print "Shake DATA  = ",$SHAKE_DAT
	print "Novak DATA  = ",$NOVAK_DAT
	print "Admit DATA  = ",$ADMIT_DAT
	print "Namaz DATA  = ",$NAMAZ_DAT
	pause

	print "解析定義のExcelから地盤情報を抽出して地盤物性を変数に代入する"
	cmd{
		excelcsv -proc ex2csv -excel analysis.xlsx -sheet <$SOIL> -csv JIBAN.csv
		cat JIBAN.csv
		}
	csvfile = JIBAN.csv
	for j = 1 to 7
		M = j + 1
		$G[j]  = csv_r[M]_c[2]
		$VS[j] = csv_r[M]_c[3]
		$RO[j] = csv_r[M]_c[4]
		$PO[j] = csv_r[M]_c[5]
		$T[j]  = csv_r[M]_c[6]
		$SH[j] = csv_r[M]_c[7]
	next j
	pause

	print "解析定義のExcelから入力波のCSVを抽出してShakeが読み込める書式に変換する"
	cmd{
		excelcsv -proc ex2csv -excel analysis.xlsx -sheet WAVE -csv WAVE.csv
		csv2wav.sh -c WAVE.csv -w <$INPUT_WAV> -r <$WNUM> -s <$WSTP>
		}
	pause

	print "Shakeデータを作成して解析を実行する"
	makedata shake.template $SHAKE_DAT
	print "Shakeのプログラムは 吉田望 氏 が作成された DYNEQ336 を使わせていただいています。"
	cmd{ sub.sh DYNEQ336 $SHAKE_DAT $SHAKE_OUT }
	pause

	print "アドミッタンスばねデータを作成して解析を実行する"
	makedata admit.template $ADMIT_DAT
	cmd{ sub.sh admit $ADMIT_DAT $ADMIT_OUT }
	pause
	
	print "Shake解析結果からNovakばね計算に必要な情報を取得する"
	cmd{
		shake2csv.sh -L $SHAKE_OUT -c SHAKE.csv
		cat SHAKE.csv
		}
	csvfile = SHAKE.csv
	for j = 1 to 7
		$SKG[j]  = csv_r[j]_c[1]
		$SKVS[j] = csv_r[j]_c[2]
	next j
	pause

	print "novakばねデータを作成して解析を実行する"
	makedata novak.template $NOVAK_DAT
	cmd{ sub.sh novak $NOVAK_DAT $NOVAK_OUT }
	pause
	
	print "アドミッタンスばねのKとCを取り出す(Cの計算は暫定的に1.0Hz)"
	cmd{
		admit2csv.sh -L $ADMIT_OUT -c ADMIT.csv -f 1.0
		cat ADMIT.csv
		}
	pause

	print "novakばねのKとCを取り出す(Cの計算は暫定的に1.0Hz)"
	cmd{
		novak2csv.sh -L $NOVAK_OUT -c NOVAK.csv -f 1.0
		cat NOVAK.csv
		}
	pause

	print "アドミッタンスばねとnovakばねの値を変数に代入する(Cの計算は暫定的に1.0Hz)"
	csvfile = ADMIT.csv
	$ADMK[1] = csv_r[1]_c[1]
	$ADMK[2] = csv_r[2]_c[1]
	$ADMC[1] = csv_r[1]_c[3]
	$ADMC[2] = csv_r[2]_c[3]
	print "SWAY K = ",$ADMK[1]," ROCK K = ",$ADMK[2]
	print "SWAY C = ",$ADMC[1]," ROCK C = ",$ADMC[2]
	csvfile = NOVAK.csv
	for j = 1 to 6
		$NVKK[j] = csv_r[j]_c[1]
		$NVKC[j] = csv_r[j]_c[3]
		print "NOVAK K= ",$NVKK[j]," NOVAK C= ",$NVKC[j]
	next j
	pause

	print "固有値解析用データ作成と実行(一度、1.0Hzの減衰で応答解析を実施する)"
	makedata namazu.template $NAMAZ_DAT
	cmd{ sub.sh namazu $NAMAZ_DAT $NAMAZ_OUT }
	pause

	print "解析モデルの地盤-建屋連成1次の振動数を取得する"
	cmd{
		namazu2csv.sh -L $NAMAZ_OUT -c NAMAZU.csv
		cat NAMAZU.csv
		}
	pause

	print "連成1次の減衰を取得する"
	csvfile = NAMAZU.csv
	$FREQ = csv_r[1]_c[1]
	cmd{
		admit2csv.sh -L $ADMIT_OUT -c ADMIT.csv -f $FREQ
		novak2csv.sh -L $NOVAK_OUT -c NOVAK.csv -f $FREQ
		}
	pause

	print "アドミッタンスばねとnovakばねの値を変数に代入する"
	csvfile = ADMIT.csv
	$ADMK[1] = csv_r[1]_c[1]
	$ADMK[2] = csv_r[2]_c[1]
	$ADMC[1] = csv_r[1]_c[3]
	$ADMC[2] = csv_r[2]_c[3]
	print "SWAY K = ",$ADMK[1]," ROCK K = ",$ADMK[2]
	print "SWAY C = ",$ADMC[1]," ROCK C = ",$ADMC[2]
	csvfile = NOVAK.csv
	for j = 1 to 6
		$NVKK[j] = csv_r[j]_c[1]
		$NVKC[j] = csv_r[j]_c[3]
		print "NOVAK K= ",$NVKK[j]," NOVAK C= ",$NVKC[j]
	next j
	pause

	print "応答解析を実施する"
	makedata namazu.template $NAMAZ_DAT
	cmd{ sub.sh namazu $NAMAZ_DAT $NAMAZ_OUT }
	pause
	
	print "応答スペクトル計算"
	pause
	cmd{
		wav2csv -in $INPUT_WAV -csv <$INPUT_WAV>.csv -fmt '(8f10.0)' -step $WSTP -dt 0.01
		wav2csv -in $NAMAZ_WAV -csv <$NAMAZ_WAV>.csv -fmt '(8f10.0)' -step $WSTP -dt 0.01
		duhamel -c <$INPUT_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 1 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 1 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 2 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 3 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 4 -b 300
		duhamel -c <$NAMAZ_WAV>.csv -d 0.05 -s 0.02 -e 5.0 -r 5 -b 300
		}
	pause

	print "一括後処理(関連するCSVファイルを全てTemplateのExcelに貼り付ける)"
	cmd{
		sr_post.sh <$CASE> <$WAVE> <$SOIL> <$BLDH> <$INPUT_WAV> <$CASE>_<$WAVE>_<$SOIL>_<$BLDH>.xlsx
		open <$CASE>_<$WAVE>_<$SOIL>_<$BLDH>.pdf
		}
	pause

next i

CalculiXの接触問題でマルチCPUの性能をテストする

接触解析サンプルのモデル図

こちらにある接触問題のサンプルで、CalculiXで使用しているSPOOLES(連立方程式ソルバ)の性能をテストしてみました。

CalculiXの並列計算に関する記述は以下です。
https://www.xsim.info/articles/CalculiX/ccx-doc/node3.html

1CPUで計算した結果と4CPUで計算した結果は以下のとおり。

1CPU4CPU
real0m21.924s0m13.749s
user0m19.113s0m28.569s
sys0m0.643s0m1.096s

実行時間で4CPUの方が1.6倍ほど速い結果となっています。
「user」が1CPUより4CPUの数値が大きくなっているのはCPU数分の数値を合算している為です。
「sys」は何故でしょう?
4CPUの方がDisk I/O が多いのでしょうか。

接触解析サンプルの結果図

◼️1CPU
Using up to 1 cpu(s) for the stress calculation.

Using up to 1 cpu(s) for the symmetric stiffness/mass contributions.

Factoring the system of equations using the symmetric spooles solver
Using up to 1 cpu(s) for spooles.

Using up to 1 cpu(s) for the stress calculation.
average force= 16.144929
time avg. forc= 12.096706
largest residual force= 0.000454 in node 1756 and dof 2
largest increment of disp= 1.282442e-01
largest correction to disp= 9.891491e-06 in node 609 and dof 2

convergence

the increment size exceeds the remainder of the step and is decreased to 0.000000e+00
Using up to 1 cpu(s) for the stress calculation.

Job finished

real 0m21.924s
user 0m19.113s
sys 0m0.643s

◼️4CPU
Using up to 4 cpu(s) for the stress calculation.

Using up to 4 cpu(s) for the symmetric stiffness/mass contributions.

Factoring the system of equations using the symmetric spooles solver
Using up to 4 cpu(s) for spooles.

Using up to 4 cpu(s) for the stress calculation.
average force= 16.144929 time avg. forc= 12.096706
largest residual force= 0.000454 in node 1756 and dof 2
largest increment of disp= 1.282442e-01
largest correction to disp= 9.891491e-06 in node 609 and dof 2

convergence

the increment size exceeds the remainder of the step and is decreased to 0.000000e+00
Using up to 4 cpu(s) for the stress calculation.

Job finished

real 0m13.749s
user 0m28.569s
sys 0m1.096s

pythonで作るconeモデル基礎底面ばね

こちらの文献を参考に作成してみました。
http://library.jsce.or.jp/jsce/open/00550/2010/16-01-0030.pdf

元のWolf氏の文献を見つけらず「z0=近似円錐地盤の頂点高さ」にどのような数値を設定すべきか不明です。

import math
import cmath
import scipy.special

G = 1.8E+07
VS = 100.0
poisson = 0.33
R = 6.0
Z0 = 15.0

for i in range(0,80):

    freq = float(i) * 0.1
    omega = 2.0 * math.pi * freq

    a0 = omega * R / VS

    K = 8.0 * G * R / ( 2.0 - poisson )
    ca0 = math.pi * ( 2.0 -poisson ) / 8.0
    a0cao = a0 * ca0
    Zsway = K * complex(1.0, a0cao) / ( G * R )

    mu = 0.3 * math.pi * ( poisson - 1.0 / 3.0 )
    Kt = 8.0 * G * R**3.0 / ( 3.0 * ( 1.0 - poisson ) )
    a02 = a0**2.0
    kta0 = 1.0 - 1.0 / 3.0 * mu / math.pi * Z0 / R * a02 - 1.0 / 3.0 * a02 / ( ( 2.0 * R / Z0 )**2.0 + a02 )
    cta0 = Z0 / ( 6.0 * R) * a02 / ( ( 2.0 * R / Z0 )**2.0 + a02 )
    Zrock = Kt * complex( kta0, a0 * cta0 ) / ( G * R**3.0  )

    print(a0, ',', Zsway.real, ',', Zsway.imag, ',', Zrock.real, ',', Zrock.imag)

pythonで作るNOVAKばね計算プログラム

novak impedance
novak impedance

建屋側面地盤のインピーダンスを簡易計算する手法としてNOVAKのばねがよく知られています。

こちらのPDFにその理論式が掲載されています。

これをもとにpythonで作成したプログラムが以下です。

import math
import cmath
import scipy.special

G = 1.0E+09
VS = 700.0
poisson = 0.3
R = 10.0
HJ = 1.0
h = 0.00

VP = VS * math.sqrt( ( 2.0 * ( 1.0 - poisson ) ) / ( 1.0 - 2.0 * poisson ) )

for i in range(1,2000):

    freq = float(i) * 0.01
    omega = 2.0 * math.pi * freq
    a0 = omega * R / VS
    b0 = omega * R / VP

    Z = cmath.sqrt( 1.0 + 2.0 * h * complex( 0.0, 1.0 ))
    a0Z = a0 * complex( 0.0, 1.0 ) / Z
    b0Z = b0 * complex( 0.0, 1.0 ) / Z
    # B = math.sqrt( ( 2.0 * ( 1.0 - poisson ) ) / ( 1.0 - 2.0 * poisson ) )
    # b0Z = a0Z / B

    # 第2種変形ベッセル関数
    k0a0 = scipy.special.kv(0, a0Z)
    k0b0 = scipy.special.kv(0, b0Z)
    k1a0 = scipy.special.kv(1, a0Z)
    k1b0 = scipy.special.kv(1, b0Z)

    p1 =  - math.pi * HJ * G * a0**2.0
    p2 = 4.0 * k1b0 * k1a0 + a0Z * k1b0 * k0a0 + b0Z * k0b0 * k1a0
    p3 = b0Z * k0b0 * k1a0 + a0Z * k1b0 * k0a0 + b0Z * a0Z * k0b0 * k0a0
    kuj = p1 * p2 / p3

    print(freq, ',', kuj.real, ',', kuj.imag)

AEiでCalculiXのBoltサンプルの荷重を変化させて歪みが一定値を超える荷重を得る

CalculiX は強力な非線形ソルバーなので、容易に静的な荷重を斬増させて非線形解析することができます。

解析コード本体で荷重を増加させていくことが可能なのですが、今回はAEi(Analysis Execution interpreter)を使って外部から静的な加力荷重を変化させて、歪みがある一定の数値を超える荷重値を調べます。

以下では、外力を0.1ずつ増加させ、歪み(Exx)が 5.0E-5 を超える荷重値を得るというデモです。

AEI and CalculiX Bolt Sample

上記は外力を線形に増加させただけですが、後処理に小さなプログラムを用意すれば、後段の解析に結果を受け渡しながら連続的に解析を実行することも可能です。

また、等価線形解析機能がないソルバーの結果を読み込んで前回解析との差分から後段解析の物性値を作成するプログラムを用意すれば、外部で等価線形解析を実施することも可能です。

MacOS(10.13.6)にCalculiXをインストールする

FreeCADでCalculiXの解析結果を表示する

CalculiXはオープンソースの有限要素法解析ソフトです。応力解析や固有値解析、過渡応答解析、そして接触問題も解くことができる強力なソルバーです。
入力ファイルの形式がAbaqusと良く似ていて、Abaqusの使用経験がある方ならデータを理解しやすいのではないでしょうか。

以下が本家サイトです。
http://www.calculix.de

WindowsやLinuxは実行形式のファイルをインストールできるインストーラが用意されていますが、MacOSはソースファイルからビルドする必要があります。
(以前はHomebrewからインストールできた様ですが現在はできません)

本家サイトにあるMacOS用のインストール資料は、Linux上でビルドする方法をMacOS用に書き換えたもので、Makefileの記述を変更したり、ソースファイルの一部を修正するように手順が示されています。

かなり面倒で、gccやfortranのバージョン、その他の環境が少し違うだけでビルドがつまづきます。

しかし、GithubにMacOS用に各種のファイルを修正してアップしてくださっている方がいます。
https://github.com/bobmel/CalculiX

このサイトからCloneをダウンロードし、本家のインストール資料に従ってビルドすればビルドが成功します。

上記のサイトでは、gcc-8とgfortran-8でコンパイルされていますが、gcc-9とgfortran-9でも問題ありませんでした。

しかし、ARPACKをビルドする際に本家の資料ではARmake.incを修正する様に記載されていますが、githubからダウンロードしたARmake.incは修正されていませんでしたので、僕は念のために以下の修正を反映しました。

LAPACKLIB = → # LAPACKLIB =
BLASLIB = → # BLASLIB =
ALIBS = $(ARPACKLIB) $(LAPACKLIB) $(BLASLIB) → ALIBS = $(ARPACKLIB)

なお、この修正をしなくてもビルドが成功するかは試していません。

C#でコマンドを実行する

以下、サンプルプログラムです。

using System;
using System.Diagnostics;

namespace test
{
    class MainClass
    {
        public static void Main(string[] args)
        {
            Process p = new Process();

            p.StartInfo.FileName = "ls"
            p.StartInfo.Arguments = "-al";
            p.StartInfo.CreateNoWindow = true;
            p.StartInfo.UseShellExecute = false;
            p.StartInfo.RedirectStandardOutput = true;
            p.StartInfo.RedirectStandardError = true;

            p.Start();
            string output = p.StandardOutput.ReadToEnd();
            string erroutput = p.StandardError.ReadToEnd();

            p.WaitForExit();
            p.Close();
            p.Dispose();

            Console.Write(output);
            Console.Write(erroutput);
        }
    }
}
MacOSでTerminal.appを起動して下記のように実行すればMacOSでも動作します。
/Library/Frameworks/Mono.framework/Versions/5.18.1/bin/mono32 test.exe