コンテンツへスキップ

画像から、全体の中心点と模様展の半径、角度を取得してファイルで出力する方法です。opencvを使用します。

画像から、全体の中心点と模様(濃い青丸)の半径、角度を取得したので、メモを修正して公開します。最終的に極座標グラフにします。

目次

  • 目的
  • 手順
  • 環境
  • コード
  • 解説

目的

最近、画像から半径と角度をたくさん取得する必要に駆られました。その作業はかなり手間なので自動化してやろうというのが目的です。

私の使用する画像をそのまま使うわけにはいかないので、似たような画像をいらすとやで拾ってきました。輪郭があんまり丸くなかったので加工しています。

螺旋の貝殻のイラスト

この画像から全体の中心点と模様(濃い青丸)の半径、角度を取得します。

手順

  1. グレースケール&位置座標を取得
  2. 円の取得の中心点の座標を取得
  3. 中心点を基準にして (rcosθ,rsinθ)で表現
  4. 各点(0~i)を(r, θ)で表記
  5. csvで出力

環境

  • Python 3.8.0
  • conda 4.8.0
  • opencv 4.1.2
  • matplotlib 3.1.2
  • numpy 1.17.3

コード

import cv2
import matplotlib.pyplot as plt
import numpy as np
import csv
import math
import sys

"""
あらかじめ同じフォルダに「white.png」と「読み込み対象.png」をおいておく必要がある。
芯の中心を原点にして各点(小円)の半径、角度を取得する。
"""
#input
input='rasen'
img = cv2.imread('input'+input+'.png',cv2.IMREAD_COLOR)#ndarray.グレイスケールで読み込みとエラー。
white = cv2.imread('white.png',0)#ndarray

#事前処理
img = cv2.medianBlur(img, 5) #メディアンフィルタを用いて画像を平滑化
cv2.imwrite("output/img.png",img)#debag1。outputフォルダに作成される。
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
cv2.imwrite("output/gray.png",gray)#debag2
edge = cv2.Canny(gray, 20, 40)
cv2.imwrite("output/edge.png",edge)#debag3
#中心円の取得
circles = cv2.HoughCircles(edge, cv2.HOUGH_GRADIENT, 1, 1200, param1=60, param2=20, minRadius = 800, maxRadius =1000)
#print(circles)
circles = np.uint16(np.around(circles))
x_c=np.int64(circles[0][0][0])
y_c=-np.int64(circles[0][0][1])
r_c=np.int64(circles[0][0][2])
center =[x_c,y_c]
norm_c= np.linalg.norm(center, ord=2)
for i in circles[0,:]:
    # 円を描く
    cv2.circle(img,(i[0],i[1]),i[2],(0,0,0),4)
#小円の取得
circles = cv2.HoughCircles(edge, cv2.HOUGH_GRADIENT, 1, 200,param1=60,param2=20, minRadius = 15 ,maxRadius = 50) #param1は色の違い、param2は形の雑さ
circles = np.uint16(np.around(circles))#整数になり、16bitに変換された配列になる。
j=0
for i in circles[0,:]:
    # 小円を描く
    j+=1
    cv2.circle(img,(i[0],i[1]),i[2],(0,0,0),4)
print(str(j)+"points")

#データの整形
#radiusは不要
radius=circles[:,:,2]
radius=np.array(radius[0,:])
x=circles[:,:,0]
x=x[0,:]
y=circles[:,:,1]
y=y[0,:]

#csvファイルの作成
xs=[]
ys=[]
points=[]
with open('output/'+input+'_x-y.csv', 'w') as f:
    writer = csv.writer(f, lineterminator='\n')
    writer.writerow(center)
    for i in range(len(x)):
        xs.append(np.int64(x[i]))
        ys.append(-np.int64(y[i]))
        points.append([xs[i],ys[i]])# points=[x,y]
        writer.writerow([xs[i],ys[i]])

# L2ノルム(距離)の取得とCSVファイルの生成
ds = []
Theta=[]
cos=[]
with open('output/'+input+'_d-Theta.csv', 'w') as f:
    writer = csv.writer(f, lineterminator='\n')
    writer.writerow(['Radius','Deg[deg]'])#ラベル
    for i in range(len(xs)):
        X = np.array([xs[i]-x_c,ys[i]-y_c])#中心からのベクトル
        norm = np.linalg.norm(X, ord=2)#中心からの距離
        ds.append(norm)
        cos.append(X[0]/norm)#cos = X[0]/SQRT(X[0]^2+X[1]^2)
        if ys[i]>y_c:#芯の中心より上側。
            Theta.append(360-math.degrees(math.acos(cos[i])))
        else:        #芯の中心より下側。
            Theta.append(math.degrees(math.acos(cos[i])))
        writer.writerow([ds[i],Theta[i]])

cv2.imwrite('output/'+input+'circle.png',img)# debag4

なお、角度は時計回りに作っています。

メモ:

画像は、色情報の2次元マトリックスになっており、画像処理の基本は左上から右上まで走査したあと一個下の行を走査するイメージです。取得した各点の順番は左上が先だと思います。

解説

1. グレースケール&xy座標で画像を取得

読み込み

cv2.imread(img,flag)
  • - cv2.IMREAD_COLOR : カラー画像として読み込む.画像の透明度は無視される.デフォルト値
  • - cv2.IMREAD_GRAYSCALE : グレースケール画像として読み込む
  • - cv2.IMREAD_UNCHANGED : アルファチャンネルも含めた画像として読み込む 上記のフラグを使う代わりに,単に1, 0, -1 の整数値を与えて指定することもできます.

画像を扱う -OpenCV より

フィルタリング

ここではあまり働いてませんが、やると精度が上がります。下の関数はC++(のはず)の記法ですが、引数の参考まで。

medianBlur(const Mat& src, Mat& dst, int ksize)

画像フィルタリング -openCV

グレイスケール

cv2.cvtColor(img, code)
  • code – 色空間の変換コード.説明を参照してください

色空間の変換 -OpenCV

debag2: グレイスケールまでの画像

エッジ処理

輪郭を取得します。

cv2.Canny

Canny -openCV

debag3 :エッジ検出までの画像

2. ハフ変換で円を検出する

ハフ変換は画像中の直線や円などを検出する操作です。ここでもpython用のサイトより見やすかったのでC++(のはず)へのリンクを張っておきます。

HoughCircles(Mat& image, vector<Vec3f>& circles, int method, double dp, double minDist, double param1=100, double param2=100, int minRadius=0, int maxRadius=0)¶

検出された円の (中心の x 座標, 中心の y 座標, 半径) のタプルを返します。

  • - image – 8ビット,シングルチャンネル,グレースケールの入力画像.
  • - circles – 検出された円を出力するベクトル.各ベクトルは,3要素の浮動小数点型ベクトル (x, y, radius) としてエンコードされます.
  • - method – 現在のところ, CV_HOUGH_GRADIENT メソッドのみが実装されています.基本的には 2段階ハフ変換 で,これについては Yuen90 で述べられています.
  • - dp – 画像分解能に対する投票分解能の比率の逆数.例えば, dp=1 の場合は,投票空間は入力画像と同じ分解能をもちます.また dp=2 の場合は,投票空間の幅と高さは半分になります.
  • - minDist – 検出される円の中心同士の最小距離.このパラメータが小さすぎると,正しい円の周辺に別の円が複数誤って検出されることになります.逆に大きすぎると,検出できない円がでてくる可能性があります.
  • - param1 – 手法依存の 1 番目のパラメータ. CV_HOUGH_GRADIENT の場合は, Canny() エッジ検出器に渡される2つの閾値の内,大きい方の閾値を表します(小さい閾値は,この値の半分になります).
  • - param2 – 手法依存の 2 番目のパラメータ. CV_HOUGH_GRADIENT の場合は,円の中心を検出する際の投票数の閾値を表します.これが小さくなるほど,より多くの誤検出が起こる可能性があります.より多くの投票を獲得した円が,最初に出力されます.
  • - minRadius – 円の半径の最小値.
  • - maxRadius – 円の半径の最大値.

cv::HoughCircles¶ -OpenCV より

出力は中心位置と半径です。単位は画素だと思われます。

メモ:

上記コードでは以下のように使用しています。

circles = cv2.HoughCircles(gray, cv2.HOUGH_GRADIENT, 1, 180, param1=60, param2=40, minRadius = 800, maxRadius =1100)

上のコードでは入力画像は2504*2607です.minRadius = 800なので、画像の半分以上を覆うような大きい円だけを検出しています。 param2は経験上、円の雑さに関するパラメータで、小さいほど歪な円を検出します。

円を描画

void circle(Mat& img, Point center, int radius, const Scalar& color, int thickness=1, int lineType=8, int shift=0)¶
  • - img – 円を描画する画像.
  • - center – 円の中心座標.
  • - radius – 円の半径.
  • - color – 円の色.
  • - thickness – 円の枠線の太さ.負の値の場合,円が塗りつぶされます.
  • - lineType – 円の枠線の種類, Line() の説明を参照してください.
  • - shift – 中心点の座標と半径の値において,小数点以下の桁を表すビット数

cv::circle -OpenCV より

debag4: ハフ変換までの画像

3. 円の中心点の座標を取得

各点は(x,y,r)を取得します。芯の位置ベクトル($x_c,y_c$)から、中心点からの各点の位置ベクトル(x_i,y_i)=(x-x_c,y-y_c)を求めます。x軸単位ベクトルを(x_e,y_e)=(1,0)とすると、内積から、

$$cos\theta = \frac{(x-x_c)x_e+(y-y_c)y_e}{\sqrt{x_c^2+y_c^2}\sqrt{(x-x_c)^2+(y-y_c)^2}}$$

$$cos\theta = \frac{x-x_c}{\sqrt{(x-x_c)^2+(y-y_c)^2}}$$

$$\theta = deg(arccos(\frac{x-x_c}{\sqrt{(x-x_c)^2+(y-y_c)^2}}))[deg]$$

【NumPy入門】ベクトルの大きさ(ノルム)を計算するnp.linalg.norm

np.linalg.norm(X, ord=2)

L2ノルム(ユークリッドノルム)は一般的な"長さ"です。

ノルムの意味とL1,L2,L∞ノルム -高校数学の美しい物語

グラフにして確認

結果が正しいかを確認します。アバウトでOKなので、画像とグラフを重ねていきます。

細かいことはかつて書きました=>[python]データのインポートと極座標系のグラフを表示- イノマタの趣味部屋

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

"""
CSVファイルから極座標系のグラフをプロットする
入力:(CSVファイル名, 出力画像ファイル名,半径,角度)
出力:出力画像ファイル名.png
"""
def polar(in_file,out_file,param_r,param_theta):    
    df = pd.read_csv("output/"+in_file)#ファイル読み込み
    r = df[param_r]
    deg = df[param_theta]
    radian = np.deg2rad(deg)#
    Len=len(r)

    fig = plt.figure(figsize=(10, 10))#新しい図面
    ax = fig.add_subplot(111, projection='polar')#軸の表示
    ax.set_theta_direction(-1)#時計回り
    #ax.plot(radian,r, '-o')
    ax.plot(radian,r, 'or')
    
    i=2
    #while i < Len:
    #    ax.plot([radian[i],radian[i-1]],[r[i], r[i-1]], 'xr-')
    #    i +=2
    plt.savefig(out_file)
    plt.show()

polar('0%_d-Theta.csv','0%.png','Radius','Deg[deg]')
取得した数値からグラフを作成

なお、画像は下記サイトで背景を透明にしました。

画像に透明色を設定する -peko step

はじめの画像と重ね合わせると、うまくグラフにできていることが確認できます。 取得したCSVは正確みたいですね。

グラフと元々の画像を重ねた画像

もともとの画像を修正したりエッジを強調したりすれはもっといい結果が得られるはずです。

雑記

パラメータの変更も自動化できれば良かったですが、できませんでした。

目的は達成したものの、手で測定するのと同じくらい時間がかかってしまったのも悲しいです。追加があれば効果が発揮できるんですが。

参考

画像処理の基本は、古い本ですがこの本が良いと思います。

コンピュータ画像処理 田村 秀行 (著) オーム社

他に、最新のアルゴリズムとか実用的な部分はOpenCVを見るのが早いでしょう。

研究などでpythonを使って計算したりグラフ描画したりしたい人のための基礎事項まとめ本」でした。

「研究などでpythonを使って計算したりグラフ描画したりしたい人のための基礎事項まとめ本」でした。

科学技術計算のためのPython入門――開発基礎,必須ライブラリ,高速化 技術評論社 中久喜 健司 (著)

期間: 2019年9月4週目~10月1週目くらい

内容(個人的まとめ)

  1. pythonプログラミングの例
  2. IPython,Spyderなど
  3. pythonの基礎事項
  4. NumPy
  5. pandas
  6. SciPy
  7. Matplotlib
  8. 高速化

特徴

科学技術計算でPythonを使用したいと考えている方に向けた本になっていました。特定の何かをしようというものではなく、例が豊富で実行しながらpythonでできることを一通り網羅的に紹介して貰える本です。

感想

(参考として)読む前:pythonのフロー制御とモジュールを使って簡単な計算を自動化するプログラムが作成できる程度の知識がありました。一方で、Open CVなどその時々に必要なモジュールがないかググり、Qiita等に書いてある範囲内で使用することしかできませんでした。

全体

読む以前から多少の知識がある状態で、基本事項を改めて確認し直したい私にはとてもちょうどいい本でした。もしかしたらpythonを書いたことのない方には難しいかもしれませんが、例が豊富なのでひとまず読むことはできるかと思います。

公式リファレンスを読むのはまだ重い...という私と同じような方におすすめです。

1. pythonプログラミングの例

プログラムを書き、テストし、高速化して利用するという流れを一通り見せてくれようとしている章でした。作業フローを教えてくれる本はあまりであったことがなかったので貴重です。

例としてロケットシミュレーターを挙げており、公開されているコードを落としてこれるので実行しながら読み進められるはずです。が!地図をプロットするmatplotlib.basemapがエラー吐きまくりだったのでコードの実行はあきらめました。

2. IPython,Spyderなど

私はJupyter Notebookで普段書いています。Jupyter Notebookで書くとグラフや画像がその場で表示されるしTab補完もできるのでで重宝しています。

IPythonとはどんなものか、Spyderとは何かといった、知らなくてもどうにかなるが知っていると便利、でもわざわざ調べないことが書いてあり有難かったです。

3. pythonの基礎事項

一応知っていた内容が多い部分でした。「インデキシングとスライシングっていう名前の上での分類があったんだー」というような、教科書改めて読むと面白いという感じの内容です。

4. NumPy

NumPyいつもimportするけどnp.piしか使わないじゃん、と思っていた俺を殴りたい。多次元配列ndarrayとはどのようなものか、何がいいのかということを、例を実行しながらイメージできました。

5. pandas

よく聞くpandasとNumPyの違いとは何か、どんな時にpandasを使えばいいのか、なんとなく分かった気になれます。主にデータ型Series,DataFrame,Panelsに関して書かれています。

ただし、Panelsに関しては、「将来的になくなるよ、使わない方がいいよ」と警告してもらったので(下記)読み飛ばしました。本の悪い部分は情報の更新がわからないことなので仕方がないですね。

Panel is deprecated and will be removed in a future version.
The recommended way to represent these types of 3-dimensional data are with a MultiIndex on a DataFrame, via the Panel.to_frame() method
Alternatively, you can use the xarray package http://xarray.pydata.org/en/stable/.

6. SciPy

SciPyの項目では、以下の項目について例を挙げていました。

  • 統計関数
  • 離散フーリエ解析
  • ボーデ線図
  • データの内挿(グラフに近似線を追加するなど)
  • デジタル信号フィルタの設計
  • 行列の分解

7. Matplotlib

使いたいときにググればいいと思い、読み飛ばしました。

8. 高速化

自分が今必要としていない部分だと思いましたので、読み飛ばしました。

次に読みたい

関連記事

pythonでグラフを表示するまでの流れを紹介します。グラフは自由振動、減衰振動、強制振動、倒立振り子(倒立振子、縦棒制御)の制御っぽいことをして時間変化をグラフで表示しました。おまけでCR回路の電圧の充電具合をプロットしました。 グラフを表示するのはCとかでも簡単にできますが、機械学習やってみたいのでpythonに触れる機会にします。

目次

  • 目的
  • グラフの様子
  • pythonとは
  • 環境とインストール
  • Jupyter Notebookの使い方
  • 振動の基礎とコード
  • 倒立振り子の運動方程式
  • 雑記
  • 参考

目的

思いのままのグラフを表示します。多項式だけではなく行列やベクトル、微分方程式も扱えるようになりたいと思います。

グラフの様子

今回は2つの似たような数式モデルをプロットしました。

図1:倒立振り子
  • 倒立振り子(縦棒制御)

手の上で箒のバランスを取るやつです。振動の中でも強制振動です。ニュートンの運動方程式を立てると、次のようになります。振り子は2l[m]、棒はm[kg]、台車はM[kg]です。このとき、角度Θを十分に小さいとして 近似してかなり簡単な式にまとめてあります。導出過程はページ最後に書きました。

$$\dot \theta_1(t)=\theta_2 \\ \dot \theta_2(t) = \frac{3}{l(4M+m)}(-L\theta_2+{g(M+m)-K}\theta_1)$$

簡単のためにl=1[m],M=2/3[kg],m=1/3[kg]を選ぶと、

$$\dot \theta_1(t)=\theta_2 \\ \dot \theta_2(t) = -L\theta_2+(g-K)\theta_1$$

この式で常微分方程式を解くと下記のような美しいグラフに!!コードが汚いのは御愛嬌ということで。python初めてだったのでほぼパクリです。減衰比については振動の基礎の項目に軽く書いています。

図2:倒立振り子の角度制御
from sympy import *      # sympyライブラリから全ての機能をimoprt
import numpy as np       #ベクトルや行列を表すのに使うライブラリnumpy
from scipy.integrate import odeint #scipyに含まれる常微分方程式を解く関数
import matplotlib.pyplot as plt #グラフを書くためのライブラリ
%matplotlib inline     #jupyter notebook内で表示させるため

"""
2階微分方程式 倒立振子 
"""
#シンボル定義
x=Symbol('x')                  # 文字'x'を変数xとして定義
t=Symbol('t')                  # 文字 't'を変数tとして定義

g = 9.8     #重力加速度
La = 0.89   #減衰比0.7のグラフのフィードバックゲイン
Ka = 10.2  
Lb = 1.0    #減衰比1.0のグラフのフィードバックゲイン
Kb = 10.05  
Lc = 0.0    #減衰比0.0のグラフのフィードバックゲイン
Kc = 10.0   
Ld = 1.0    #減衰比0.3のグラフのフィードバックゲイン
Kd = 12.57  

def func1(x, t, K):
    x1,x2=x
    dxdt=[x2,(-La*x[1])+((g-Ka)*x[0])]
    return dxdt
#同じことを複数書いているので省略

x0 = [10.0,2.0] # 初期条件 [x(0), dx(0)/dt]を表す

t=np.linspace(0,30,1000) # 時間 0から10までを101等分する
sol1=odeint(func1, x0, t, args=(k,)) # 微分方程式を解く 
#同じことを複数書いているので省略

#グラフ表示
plt.plot(t, sol3[:,0], linewidth=1,label='ζ=0.0') #  角度を図示。
plt.xlabel('time t[s]', fontsize=18)
plt.ylabel('angle θ[rad]', fontsize=18)
plt.legend(loc='upper right')
plt.savefig('huriko1.jpg')
plt.show()
  • LCR回路

そしてもう一つの式がこれです。回路のキャパシタはダンパー、インダクタはばねと同じです。

$$ E = v_R + v_C + v_L \\ E = L\frac{di(t)}{dt}+Ri(t)+\frac{1}{C}\int i(t)dt $$

これをキャパシタ電圧でまとめると

$$ i = C\frac{d}{dt}v_C \\ \frac{d^2}{dt^t} v_c = -\frac{R}{C} \frac{d}{dt}v_c -\frac{1}{LC}v_c + \frac{E}{LC}$$

これは強制振動の式と同じ形なので、同じようなグラフがかけます。

図3:LRC回路のVcの時間応答

python とは

今回のpythonは読みやすく、それでいて効率もよいコードをなるべく簡単に書けるように作られたプログラミング言語です。言語自体も書きやすいように作られていますが、ライブラリやワーススペースが豊富で、だいたい何でも出来るみたいです。最近ではAIや機械学習が人気ですが、そういう事をするためのライブラリが豊そろっているみたいで、最近活躍している言語です。 2系と3系は全然違うので注意が必要です。

ライブラリというのは、

ライブラリとは、複数のプログラムが共通して利用するコードをまとめたファイルである。 ライブラリ自体は単独で実行することはできない。 from I-8-8. ライブラリとプログラムの関係 | 日本OSS推進フォーラム

というものです。必要なときに呼び出すと自分の仕事だけやって帰って行きます。

ちなみに、初学者へ。世の中のホームページやアプリの画面表示は別な言語で書いています。pythonとかCとかは計算処理をやってくれるものなので、これだけやっても画面表示は出来ません。注意してください。勉強しはじめたころの私はこのあたりよく分かって無かったために時間が掛かってしまったので一応書いておきました。

環境とインストール

環境はAnacondaでまとめてインストールしてJupyter Notebookで書きます。 python公式でインストールしても良いんですが、他にもいろいろまとめてインストール出来るらしいDownloads - Anacondaでインストールした方が楽みたいです。デフォルトのままで良いので指示に従うだけです。

Jupyter Notebookの使い方

インストールしたらJupyterNotebookを起動します。するとコマンドコマンドプロンプト(っぽいなにか)が出てきますので、そこに書かれているURLをブラウザで開いてください。

ユーザーのフォルダの中身が表示されますので、好きなのを選んで右上のnewからpython3をクリック。すると、プログラムを書くためのセルが出てきます。

今数字が書いてある色塗りの部分を押すとコマンドモードになり数式やらがかけます。In[]やOut[]と書いている部分はエディットモードで、コピーやらが出来ます。エディットモードでHを押すとヘルプが出ます。

コマンドモードで上のcodeの部分をmorkdownにすると、morkdownで文章を書けます。 morkdownは私がこのブログの下書きやメモに普段使いしているマークアップ言語(HTMLの仲間)で、# ABC と書くとタイトル的な大きな文字、## ABCと書くと1段階小さい文字、なにもつけないと普通の文字、- と書くと箇条書き、というようにかなり使いやすい言語です。

pythonの記法はすっ飛ばします。numpyというライブラリを読みこんでベクトル(1次元配列)や行列(2次元配列)を表せます。 スライシング。bool。

コマンドプロンプトでサーバーをlocalhostで立てて、ブラウザで表示しているんじゃないかと思います(このあたりよく解ってませんが。たぶん。。なので、ブラウザで開いている間にJupyterNotebookを消すとエラーをはきます。

発生したエラーについて

  • 保存した画像が真っ白になっている。jupyterにはグラフが出力されている

Python matplotlibのグラフが保存できない(windows+anaconda) teratail

振動の基礎とコード

図4:線形振動系の力学モデル

機械屋さんにとってのかなり基本的なモデル(らしい)で考えます。ニュートンの運動方程式を立てます。

$$ f(t)= m\frac{d^2x(t)}{dt^2}+c\frac{dx(t)}{dt}+kx(t) $$

これは強制振動を表しています。

振動は、自由振動、減衰振動、強制振動の3つに分けられます。自由振動は正弦波などのように振動し続けます。減衰振動は自由振動が次第に減衰していきます。強制振動は減衰振動に外力を加えている場合です。 日常生活では大概空気抵抗などで減衰するので減衰振動がイメージしやすいかと思います。一方、制御するときはモータが動いているので強制振動です。

外力をf(t)と書いていますから、f(t)=0なら減衰振動、f(t)の項が有れば強制振動になります。

上の式を作為的に以下のように変形します。

$$ f(t)= \frac{d^2x(t)}{dt^2}+2\zeta\omega_n\frac{dx(t)}{dt}+\omega_n^2 \ \\減衰比\zeta = \frac{c}{2\sqrt{km}}, 固有振動数\omega_n = \sqrt{\frac{k}{m}} $$

減衰比ζは減衰具合、固有振動数ω_nは振動時の周波数に対応しています。これが減衰振動の基本的なモデルです。ζ=0のとき全く減衰しないのでsin波のような自由振動をし、ζ=1の時過不足無く(クッションに沈むように)目標値に到達します。実際には誤差や外乱が入るので制御系を作るときには0.7くらいにするそうです。

今回はpythonで解きたいので、ここから変形していきます。本当はいろんな方法が有るみたいなんですが、そのうちやると言うことにして、今回は参考本と同じ方法で解いていきたいと思います。

常微分方程式を解きます。二次方程式は大変なので、以下のように書きかえます。

$$\frac{dx_1(t)}{dt}=x_2 \\ \frac{dx_2(t)}{dt} = -\frac{c}{m}x_2-\frac{k}{m}x_1 +f(t)$$

これを行列を用いたベクトル表現にすると、

$$ \frac{d}{dt} \left[ \begin{array}{rrr} x_1 \\ x_2 \end{array} \right] = \left[ \begin{array}{rrr} 0 \ & 1 \ \\ -\frac{k}{m} & -\frac{c}{m} \\ \end{array} \right] \left[ \begin{array}{rrr} x_1 \\ x_2 \end{array} \right] + \left[ \begin{array}{rrr} 0 \\ 1 \end{array} \right] f(t)$$

この式のように変数の位置x(t)を配列としてpythonで定義すると、以下のようになります。

m = 1.0 #質量
c = 0.5 #粘性減衰
k = 1.0 #ばね定数

def func(x, t, k):
    x1,x2=x
    dxdt=[x2,-c/m*x[1]-k/m*x[0]+sin(t/2)]
    return dxdt

x0 = [0.0,0.0] # 初期条件 [x(0), dx(0)/dt]を表す.今は初速0.

t=np.linspace(0,50,100000) # 時間 0から10までを101等分する
sol=odeint(func, x0, t, args=(k,)) # 数値的に微分方程式を解き, x(t)とx'(t)をsolのリストの[:,0]、[:,1]成分へ格納する。

以上のような感じで、パラメータをいじることで自由振動、減衰振動、強制振動がそれぞれ再現出来ました。

  • 自由振動(図5) 初期x=0.0, 初期x' = 1.0, m = 1.0, c = 0.0, k = 1.0
  • 減衰振動(図6) 初期x=0.0, 初期x' = 1.0, m = 1.0, c = 0.5, k = 1.0
  • 強制振動(図7) 初期x=0.0, 初期x' = 0.0, m = 0.0, c = 0.5, k = 1.0 f(t) = sin(t/2)
  • 強制振動(図8) 初期x=0.0, 初期x' = 0.0, m = 1.0, c = 0.5, k = 1.0 f(t) = 1
図5:自由振動
図6:減衰振動
図7:強制振動
図8:強制振動

倒立振り子の運動方程式

再度掲載

台車と振り子を図に示します。 粘性摩擦係数をB、 軸の粘性摩擦係数をCとおきました。ここから、下の4つの式が立てられます。

振り子について

$$J\ddot \theta = Vlsin\theta - Hlcos \theta - C\dot \theta …(1)\\ (J=\frac{ml^2}{3})\\ m\frac{d^2}{dt^2}(x + lsin\theta)= H…(2) \\ m\frac{d^2}{dt^2}(lcos\theta)=V - mg…(3)$$

台車について

$$M \ddot x = f(t) - B \dot x - H…(4)$$

簡単のためにC = 0, B = 0とします。ここからH,Vを消去します。

$$H = m \ddot x + ml \ddot \theta sin\theta - ml (\dot \theta)^2 sin\theta…(5) \\ V = -ml \ddot \theta sin \theta - ml(\dot \theta) cos \theta +mg…(6)$$

これを(1)(4)に代入してまとめると、

$$ J\ddot \theta= -ml^2 \ddot \theta - lmcos\theta \ddot x + mglsin\theta…(7) \\ (M+m) \ddot x = -ml\ddot \theta cos\theta + ml(\dot \theta)^2 sin\theta+f(t)…(8)$$

この式からxの2階微分を消去して、Θが十分小さいとして線形化すると

$$cos\theta\approx 1,sin\theta\approx\theta, (\theta)^2 \approx 0 \\ \ddot \theta l (4M+m)-3(M+m)g \theta = -3f(t)…(9)$$

f(t)に角度と各速度についてフィードバックさせます。

$$ f(t)= L \dot \theta + K \theta\ $$

簡単のために l = 1[m], M = 2/3[kg], m = 1/3[kg]とおくと、以下のようにかけます。プログラムで解けるように一階化します。

$$ \dot \theta_1 = \theta_2 \\ \dot \theta_2 = -L \theta_2 + (g-K)\theta_1$$

一応全体の全体をラプラス変換して伝達関数の式を求めます。(せっかく導出したから書きたい。)

$$L{\theta(t)}=\Theta(s),L{f(t)}=F(s)として\\ \frac{\Theta(s)}{F(s)}=\frac{-3}{(4M+m)ls^2+3ls-3(M+m)g+K}$$

雑記

プログラミング系のサイトは最近公開されたものを見るように気をつけるようになりました。機械系や電気系の、特に理論はかなり成熟しているのでふるいサイトでも何でも読みやすくわかりやすければ良いんですが、情報系はツールがどんどん更新されていくので最新のを見ないと妙なところで躓いてしまいます。 おかげで私のPCちゃんはパスとかアプリとか訳が分からない感じになっています。情報系の特にweb系のブロガーさんは情報収集がすごいので、最新のでも誰かしら先にやっていますので安心ですね。先人は偉大です。 それと、解析力学が難しいですね。もう少し時間をかけないとわからないです。ラグランジュ難しい。制御はやって行きたいので、のんびり本を読んでいきます。
追記(2019/3/20)参考リンクをリストにしました。Amazonリンクを増やしました。

参考