手書きモデルで固有値解析

以前、試した Google の機械学習サービス Google Cloud Vision API と iPhone の Pythonista3 を使って手書きモデルから固有値解析を実行するアプリを作成してみました。

動作の雰囲気は以下です。

Python ソースは以下のとおり。

# -*- coding: utf-8 -*-
import appex
import numpy as np
from matplotlib import pyplot as plt
from requests_toolbelt.multipart.encoder import MultipartEncoder

import base64
import json
from requests import Request, Session
from io import BytesIO
from PIL import Image

#PILで開いた画像をbase64形式に変換する
def pil_image_to_base64(pil_image):
    buffered = BytesIO()
    pil_image.save(buffered, format="PNG")
    str_encode_file = base64.b64encode(buffered.getvalue()).decode("utf-8")
    return str_encode_file


def pil_image_to_base642(fp):
    str_encode_file = MultipartEncoder(fields={'encoded_image': (fp.name, fp)})
    return str_encode_file


#PILで開いた画像をCloud Vision APIに投げる
def recognize_image(pil_image):
        str_encode_file = pil_image_to_base64(pil_image)
        str_url = "https://vision.googleapis.com/v1/images:annotate?key="
        str_api_key = "取得したAPIキー"
        str_headers = {'Content-Type': 'application/json'}
        str_json_data = {
            'requests': [
                {
                    'image': {
                        'content': str_encode_file
                    },
                    'features': [
                        {
                            #'type': "TEXT_DETECTION",
                            'type': "DOCUMENT_TEXT_DETECTION",
                            'maxResults': 10
                        }
                    ]
                }
            ]
        }

        obj_session = Session()
        obj_request = Request("POST",
                              str_url + str_api_key,
                              data=json.dumps(str_json_data),
                              headers=str_headers
                              )
        obj_prepped = obj_session.prepare_request(obj_request)
        obj_response = obj_session.send(obj_prepped,
                                        verify=True,
                                        timeout=60
                                        )

        if obj_response.status_code == 200:
            with open('data.json', 'w') as outfile:
                json.dump(obj_response.json(), outfile)
                text = get_fullTextAnnotation(obj_response.text)
            return text

        else:
            return "error"
#返ってきたjsonデータの"fullTextAnnotation"部分のテキストを抽出する
def get_fullTextAnnotation(json_data):
    text_dict = json.loads(json_data)
    try:
        text = text_dict["responses"][0]["fullTextAnnotation"]["text"]
        return text
    except:
        print(None)
        return None

# 手書き画像を共有Extensionで取得する
pil_image = appex.get_image()

# Google Cloud Vision API に画像を投げて取得したデータから重量とバネ剛性を取得する
res = recognize_image(pil_image)
sp = res.splitlines()

m = list()
k = list()

print('◼️Google Cloud Vision API の認識結果')
for v in sp:
    v = v.strip()
    v = v.replace(' ', '')
    v = v.replace(',', '.')
    d = v.split('=')
    print(v, d)
    if d[0] == 'm' or d[0] == 'M':
        m.append(float(d[1]))
    if d[0] == 'k' or d[0] == 'K':
        k.append(float(d[1]))

print("◼️バネマス情報")
for v in m:
    print('m =',v)
for v in k:
    print('k =',v)

M1 = len(m)
M2 = M1 + 1
M = np.zeros((M2, M2))
K = np.zeros((M2, M2))

# 質量マトリクスと剛性マトリクスを作成
for i in range(len(m)):
    M[i][i] = m[i] / 9.80665

for i in range(len(k)):
    i1 = i
    i2 = i + 1
    K[i1][i1] = K[i1][i1] + k[i]
    K[i2][i2] = K[i2][i2] + k[i]
    K[i1][i2] = K[i1][i2] - k[i]
    K[i2][i1] = K[i2][i1] - k[i]

# 最後のばねの端部が固定点として行と列を削除して拘束条件とする
M = np.delete(M, M1, 0)
M = np.delete(M, M1, 1)
K = np.delete(K, M1, 0)
K = np.delete(K, M1, 1)

# 質量マトリクスの逆行列を計算
M_inv = np.linalg.inv(M)

# 固有値と固有ベクトルを計算
omega, v = np.linalg.eig(np.dot(M_inv, K))

# 固有値の順番を昇順にソート
freq = list()
omega_sort = np.sort(omega)
for o in omega_sort:
    freq.append(np.sqrt(o)/(2.0*np.pi))

# 固有値のソート時のインデックスを取得
# ⇒固有ベクトルと対応させるため
sort_index = np.argsort(omega)

# 固有値に対応する固有ベクトルをソート
v_sort = []
for i in range(len(sort_index)):
    v_sort.append(v.T[sort_index[i]])
v_sort = np.array(v_sort)

print('◼️固有値解析結果')
# 結果をコンソールに表示
for i in range(len(omega_sort)):
    s = '{0:d}次 {1:10.3f} Hz'.format(i+1,freq[i])
    print(s)
    for j in range(0,v_sort.shape[1]):
        s = '{:10.3f}'.format(v_sort[i][j])
        print(s)

# グラフ化のために自由度軸を作成
dof = np.linspace(len(sort_index), 0, len(sort_index)+1)

# ここからグラフ描画
# フォントの種類とサイズを設定する。
plt.rcParams['font.size'] = 9
plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。
fig = plt.figure(figsize=(3,6))
ax1 = fig.add_subplot(111)
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。
ax1.set_xlabel('Eigen Vector')
ax1.set_ylabel('Degree of freedom')

# データの範囲と刻み目盛を明示する。
ax1.set_xticks(np.arange(-2, 2, 0.5))
ax1.set_yticks(np.arange(0, dof[0]+1, 1))
ax1.set_xlim(-1, 1)
ax1.set_ylim(0, dof[0])

# データプロット 固有ベクトルの形を見る
for i in range(len(sort_index)):
    eigen_vector = np.concatenate([v_sort[i], [0]])
    s = 'Mode' + str(i+1) + ' {:10.3f}Hz'.format(freq[i])
    ax1.plot(eigen_vector, dof, lw=1, marker='o', label=s)

fig.tight_layout()
plt.legend(bbox_to_anchor=(0, 0.1), loc='lower left')

# グラフを表示する。
plt.show()
plt.close()

固有値解析のソースはWATLABさんのこちらを参考にさせていただきました。

応答スペクトルの減衰値の違いで、加速度応答の大小関係が逆転する現象

加速度応答スペクトル図

良く見かける基盤波の加速度応答スペクトルの図です。詳しく見ると、周期0.15秒の加速度応答が、減衰1%と3%で値の大小関係が逆転しています。

一般に入力が同じなら構造物の減衰が1%と3%とでは、1%の方が応答は大きくなります。

しかし、上記の様なスペクトルの特徴を持つ基盤波を入力とした場合、構造物の卓越周期が0.15秒ならば、減衰が3%の方が1%の場合より応答が大きくなることがあります。

構造物の応答は全てのモードの重ね合せであるため、周期0.15秒のモードだけで応答が決まることはありませんが、鉄塔の様に単純な構造物の上下動解析では、1次モードだけが励起されることもあり、この周期が上記の入力の0.15秒に合致してしまった場合は減衰3%の方が応答が大きくなるという不思議な現象が起こります。

長く振動解析をしていますが、過去にこの様な現象に一度だけ出くわしたことがあります。

※上記の地震波は以下のサイトからダウンロードしたものを少々編集したものです。
https://www.bcj.or.jp/download/wave/

二つの世界を渡り歩く

振動解析に関する持論です。

時刻歴データに係数を乗じて足し合わせるプログラムを作れば時間領域を制圧できます。

伝達関数データに係数を乗じて乗算除算するプログラムを作れば周波数領域を制圧できます。

そして、ホワイトノイズに精通すれば二つの世界を渡り歩く事ができます。

動的解析コードで静的解析

動的解析コードは通常、静的解析も実行できますが、稀に静的解析を実行できない場合があります。

その時は、重量と減衰をゼロとすることで静的解析を実行することができます。
例えば、動的解析の解法がNewMarkβ法であったとすると、僕はこの静的解析を「NewMarkβ法静的解析」と呼んだりします…。