2023-01-01から1年間の記事一覧

第1種チェビシェフフィルタの伝達関数の分母多項式の係数の計算

clear; close all; s = tf('s'); xis = [0.1 0.5 1]; es = sqrt(10.^(xis/10)-1); e = es(1); tmp = sqrt(1+e^(-2))+e^(-1); for n = 1:10 D = 1; for k=1:n sinhv = 1/2*(tmp^(1/n)-tmp^(-1/n)); coshv = 1/2*(tmp^(1/n)+tmp^(-1/n)); sigma_k = -sin((2*k-…

チェビシェフ多項式の係数の計算

clear; close all; syms w; chebyshev = @(c1,c0) 2*w*c1 - c0; c0 = 1; c1 = w; cs = [c0 c1]; for n=2:10 cn_prev = cs(n-1); cn = cs(n); cn_next = chebyshev(cn, cn_prev); cs(end+1) = cn_next; end figure(); hold on; for n=1:4 cnw = cs(n+1); x = …

バタワースフィルタの伝達関数の分母多項式の係数の計算

clear; close all; s = tf('s'); Ns = 1:9; for idx = 1:length(Ns) n = Ns(idx); D = 1; for jdx=1:n j = jdx; sigma_j = -sin((2*j-1)/(2*n)*pi); omega_j = cos((2*j-1)/(2*n)*pi); sj = sigma_j + i*omega_j; D = D*(s-sj); end coef(D) end function n …

管用ねじの種類

Gネジ:Gas-pipe thread(機械結合用・平行・オス/メスねじ) Rネジ:Rohrgewinde(管用・テーパー・オスねじ) Rcネジ:Rohrgewinde cone(管用・テーパー・メスねじ) Rpネジ:Rohrgewinde parallel(管用・平行・メスねじ)Rohr=Pipe Gewinde=Screw/Thr…

JIS鉄鋼まとめ

SS: 一般構造用圧延鋼材(Structure) SM: 溶接構造用圧延鋼材(Marine) SB: ボイラ及び圧力容器用鋼板(Boiler) SPC: 冷間圧延鋼板および鋼帯(Plate Hot) : SPCC(Commercial, 一般用), SPCD(Deep Drawn: 絞り用), SPCE(Deep Drawn Extra: 深絞り用) S-C…

Intel Macbook ProでUSB CポートからのHDMIが出力されない

このあたりを試す developer.apple.com USBアクセサリの許可設定が疑わしい discussions.apple.com なんと問題は結局HDMIケーブルだった。 Macbook proのGPUが使っていると思われるAMD Freesyncという機能が関係してるのかとも思ったけど、違うのか?

アンセンテッドカルマンフィルタ(UKF)

clear; close all; f = @(x) [x(:,1).*cos(x(:,2)), x(:,1).*sin(x(:,2))]; F = @(x) [cos(x(2)), -x(1)*sin(x(2)); sin(x(2)), x(1)*cos(x(2))]; n = 2; mu = [10, pi/2]; sigma = [50 1; 1 0.025]; N = 500; x = mvnrnd(mu,sigma,N); y = f(x); ylin = one…

外付けメモリの.Trashesフォルダが消せない

ターミナルでsudo rm -rf .Trashes/*しても以下のエラーが出る。 rm: .Trashes: Operation not permittedターミナル(あるいはiTerm)にFull Disk Accessを与える必要があるapple.stackexchange.com

RPiカメラをPythonから操作する(Ubuntu on Raspberry Pi)

Pythonのpicameraモジュールは3年前に更新がとまっており使えなかったPythonのOpenCVモジュールから普通に動画が見れた sozorablog.com 試したコード import cv2 import numpy as np cap = cv2.VideoCapture(0) fmt = cv2.VideoWriter_fourcc('m', 'p', '4'…

有制約多目的最適化問題で妥協解を求める

t12の計算時にペンローズ擬似逆行列を用いるのがポイント。 パレートフロントはSTOMで希求水準をたくさん設定して求めるのが正攻法なのだろうが、横着してMatlabのboundary()コマンドを利用した。そのため実際はパレート最適値ではない点も拾ってグラフに描…

多目的最適化におけるSTOMを利用した妥協解の計算

clear; close all; x = -2:0.001:2; f1 = @(x) exp(-x)+1.4*exp(-x.^2); f2 = @(x) exp(x)+1.4*exp(-x.^2); figure(1); hold on; % plot(f1(x),f2(x),'o'); x1 = 0:0.01:7; y1 = interp1(f2(x),f1(x),x1,'spline'); plot(x1,y1,'-','DisplayName','Pareto Fr…

多目的最適化におけるトレードオフ比を考慮した妥協解の計算

clear; close all; x = -2:0.1:2; f1 = @(x) exp(-x)+1.4*exp(-x.^2); f2 = @(x) exp(x)+1.4*exp(-x.^2); plot(f1(x),f2(x),'o'); hold on; x1 = 0:0.1:7; y1 = interp1(f2(x),f1(x),x1,'spline'); plot(x1,y1,'-'); xlim([0,7]); ylim([0,7]); syms f1_symb…

満足化トレードオフ法(STOM)によるパレートフロントの計算

その1 clear; close all; aLevels = [2 1; 1 3; 0.4 0.5]; f1 = @(x) x.^2; f2 = @(x) (x-2).^2; x=-10:0.01:10; f1i = min(f1(x)); f2i = min(f2(x)); c1 = 2; c2 = 2; f1a = f1i + c1; f2a = f2i + c2; xopts = []; w1s = []; w2s = []; for idx=1:length…

基準点法によるパレートフロントの計算

clear; close all; % aspiration levels f1a,f2a aLevels = [1 3; 0.5 0.5; 0.1 0.8 2 1]; f1 = @(x) x.^2; f2 = @(x) (x-2).^2; c = 1; x = -10:0.1:10; xopts = []; for idx=1:length(aLevels) aLevel = aLevels(idx,:); f1a = aLevel(1); f2a = aLevel(2)…

重み付きlpノルム法によるパレートフロントの計算

その1 p=2の場合 clear; close all; p = 2; f1 = @(x) x.^2; f2 = @(x) (x-2).^2; w1s = 0:0.1:1; w2s = 0:0.1:1; x = -10:0.1:10; xopts = []; for idx=1:length(w1s) w1 = w1s(idx); for jdx=1:length(w2s) w2 = w2s(jdx); F = ((w1*f1(x)).^p+(w2*f2(x))…

イプシロン制約法でパレートフロントを計算する その2

clear; close all; global epsilon; a = 1/sqrt(2); f1 = @(x) 1-exp(-(x(:,1)-a).^2-(x(:,2)-a).^2); global f2; f2 = @(x) 1-exp(-(x(:,1)+a).^2-(x(:,2)+a).^2); x1opts = []; x2opts = []; A = []; b = []; Aeq = []; beq = []; x0 = [0,0]; lb = [-4,-4…

ε-制約法でパレートフロントを計算する

clear; close all; clear; close all; global epsilon; f1 = @(x) x.^2; global f2; f2 = @(x) (x-2).^2; xopts = []; A = []; b = []; Aeq = []; beq = []; x0 = [0]; lb = [-10,-10]; ub = [10,10]; epsilons = linspace(0,4,20); for idx = 1:length(epsi…

Matlabで制約つき非線形最適化問題を解く

fmincon()を使う jp.mathworks.com@nonlconの「@」がキモ clear; close all; f = @(x) x(1)^2+x(2)^2; A = []; b = []; Aeq = []; beq = []; x0 = [0,0]; lb = []; ub = []; fmincon(f, x0, A, b, Aeq, beq, lb, ub, @nonlcon) function [c,ceq] = nonlcon(…

線形加重和法によるパレート最適値の計算

clear; close all; x1 = -4:0.1:4; x2 = -4:0.1:4; f1 = 1-exp(-(x1-1/sqrt(2)).^2-(x2-1/sqrt(2)).^2); f2 = 1-exp(-(x1+1/sqrt(2)).^2-(x2+1/sqrt(2)).^2); plot(f1,f2,'rx'); hold on; w1s = 0:0.1:1.0; w2s = 0:0.1:1.0; x1opts = []; x2opts = []; for …

線形加重和法によるパレートフロントの計算

clear; close all; x = -10:0.1:10; f1 = x.^2; f2 = (x-2).^2; plot(f1,f2,'rx'); hold on; xlim([-0.2 0.2]); w1s = 0:0.1:1.0; w2s = 0:0.1:1.0; xopts = []; for i=1:length(w1s) w1 = w1s(i); for j=1:length(w2s) w2 = w2s(j); C = w1*f1+w2*f2; xopt …

増速比V/vと移行軌道の関係

clear; close all; r = 1; V_v = [1 1.1 1.2 1.3 1.31 1.4]; R = r./(2*(1./V_v).^2-1); es = (R-r)./(R+r); as = (R+r)/2; bs = as.*sqrt(1-es.^2); theta = linspace(0,2*pi,100); for idx=1:length(bs)-1 b = bs(idx); a = as(idx); e = es(idx); x = a*c…

楕円の方程式

clear; close all; es = [0 0.2 0.5 0.8 0.999]; a = 1; L = 2*a; for idx=1:length(es) e = es(idx); D = L*e; c = D/2; b = sqrt(a^2-c^2); theta = linspace(0,2*pi,50); x = a*cos(theta)-c; y = b*sin(theta); plot(x,y,'DisplayName',sprintf('e=%.3g'…

パラメタ励振のシミュレーション(マシュー方程式)

clear; close all; a = 6; q = 2.939; tspan = [0 100]; y0 = [0 0.01]; [t,y] = ode23s(@(t,y) odefun(t,y,a,q), tspan, y0); plot(t,y(:,1),'-r',t,y(:,2),'-b'); legend('x','v'); big; function dydt = odefun(t,y,a,q) dydt = zeros(2,1); dydt(1) = y(…

減衰強制振動の過渡応答

clear; close all; m = 1; k = 1; wn = sqrt(k/m); wd = wn; c = 0.1; F = 1; zeta = c/(2*sqrt(m*k)); beta = zeta*wn; t = linspace(0,100,500); ws = [1,0.2,5]; for idx=1:length(ws) w = ws(idx); A = 1/wd*(m*w^2-k+beta*c)*w/((-m*w^2+k^2)^2+c^2*w^2…

加振周波数ωと変位の応答倍率Md/位相角Φの関係

clear; close all; zetas = [1,0.5,0.2,0.1,0.01]; w_wn = linspace(0,2,100); Mds=[]; figure(); hold on; for idx=1:length(zetas) zeta = zetas(idx); Md = 1./sqrt((1-w_wn.^2).^2+4*zeta^2.*w_wn.^2); plot(w_wn, Md,'DisplayName',sprintf('zeta=%.2g'…

単純支持・集中荷重のはりの曲げ:荷重点と最大たわみ点は異なる

clear; close all; l = 5; a =4.5 disp('荷重作用点'); b = l-a; x = sqrt(1/3*(l^2-b^2)) disp('最大たわみ点'); dx =linspace(0,l,20)'; dy =zeros(length(dx),1); plot(dx,dy,'r-'); hold on; plot(x,0,'go'); plot(a,0,'k*'); x1 =linspace(0,a,20)'; x2…

PandasとRのデータフレーム操作の比較

pythondatascience.plavox.info

pecoで最近訪問したフォルダに移動する

C-x C-rに割り当てた .zshrc function peco-visit-last-folder() { local selected_dir=$(dirs -v | cut -f 2- | peco --query "$LBUFFER") if [ -n "$selected_dir" ]; then BUFFER="cd ${selected_dir}" zle accept-line fi zle clear-screen } zle -N pec…

時系列解析のためのブックガイド

logics-of-blue.com

Rでタグチメソッド(その2)

自分で実装してみた 得られた気づき: 1. 直交表は人のデータを再現するときは、そのデータに合う直交表を自分で作成する。 2. 自分で実験計画をデザインするときは、直交表のジェネレータを使って自分で適当に設定する。2水準のときはFrF2パッケージのFrF2(…