2023-05-01から1ヶ月間の記事一覧

満足化トレードオフ法(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'…