シミュレーションのイ​ンプットパラメータを​最小化するには？

例えば、物体の質量・初速・射出角度を与えて、重力と空気抵抗の影響を考慮した運動計算を行い、到達高度・落下位置までの距離を出力するシミュレーションにおいて、所定の条件を満たすように質量を最小化する（設計値を目的関数に据える）ことは可能でしょうか？ ＞x=質量、初速、射出角度, 非線形制約１；到達高度　設定値以上, 非線形制約２；落下位置～射点　設定値以上,　目的関数=質量　　の設定では、初期値付近で最適化を終了してしまいます。 ＞x=質量、初速、射出角度, 非線形制約１；落下位置～射点　設定値以上,　目的関数=到達高度-設定値　として、質量の初期値の与え方を工夫するのが一般的なのでしょうか？

[Jiro Doke]

質問からは少しイメージしづらいので、試されたコードがありましたら、再現できるものを追加して下さい。

Answer 1

[Tohru Kikawada]

Jiro Dokeさんのコメントのとおりもう少し具体的な処理がイメージできる形で書いていただけると迅速な回答が期待できます。

参考： Tutorial: MATLAB Answersで早く的確な回答を得るためのポイント

いただいた情報を元に簡単なサンプルを作ってみました。

一般的にこのような非線形最適化問題を fmincon のような勾配ベースの探索手法を使うと局所解に陥ってしまう可能性があります。

大域的な探索を行いたい場合には Global Optimization Toolbox の最適化ソルバーをご活用ください。

function [x,f,eflag,outpt] = ans_324101

xLast = []; % Last place computeall was called
myf = []; % Use for objective at xLast
myc = []; % Use for nonlinear inequality constraint
myceq = []; % Use for nonlinear equality constraint
x_pts = [];
y_pts = [];

fun = @objfun; % the objective function, nested below
cfun = @constr; % the constraint function, nested below

options = optimoptions('fmincon','MaxFunctionEvaluations',500,...
    'PlotFcn',@optimplotfval,'OutputFcn',@outfun);

% Solve Second-Order ODE with Initial Conditions
% syms y(t) m k y0 g theta v0;
% Dy = diff(y);
% y(t) = dsolve(m*diff(y,t,2) == -k*diff(y,t) - m*g,y(0)==y0, Dy(0)==v0*cos(theta))
% myf = matlabFunction(y,'file','my_y','vars',{t,g,v0,y0,k,m,theta})

% Initial values
x0 = [10 deg2rad(30) 1];

% Conditions
x_0 = 0;
y_0 = 0;
y_threshold = 3;
x_threshold = 3;

% Visualization
figure;
hopt = plot(0,0);
hold on;
plot([x_0+x_threshold; x_0+x_threshold],[y_0;  y_0 + y_threshold*1.5],'r');
plot([x_0; x_threshold*1.5],[y_0 + y_threshold;  y_0 + y_threshold],'r');
axis([0 x_threshold*1.5 0 y_threshold*1.5]);

% Lower and upper bounds
lb = [0 0 0];
ub = [10 deg2rad(90) 10];

% Call fmincon
[x,f,eflag,outpt] = fmincon(fun,x0,[],[],[],[],lb,ub,cfun,options);

      function y = objfun(x)
          if ~isequal(x,xLast) % Check if computation is necessary
              [myf,myc,myceq,~,x_tmp,y_tmp] = myModel(x);
              xLast = x;
              x_pts = x_tmp;
              y_pts = y_tmp;
          end
          % Now compute objective function
          y = myf;
      end

      function [c,ceq] = constr(x)
          if ~isequal(x,xLast) % Check if computation is necessary
              [myf,myc,myceq,~,x_pts,y_pts] = myModel(x);
              xLast = x;
          end
          % Now compute constraint functions
          c = myc; % In this case, the computation is trivial
          ceq = myceq;
      end

      function [f,c,ceq,t,x_pts,y_pts] = myModel(x)
          % Ref: https://ja.wikipedia.org/wiki/%E6%96%9C%E6%96%B9%E6%8A%95%E5%B0%84
          v_0 = x(1);
          theta = x(2);
          m = x(3);
          t = 0:0.01:5;
          k = 0.5;
          g = 9.8;
          x_pts = m*v_0/k*(1-exp(-k/m*t))*cos(theta)+x_0;
          y_pts = -m/k*((v_0*sin(theta) + m/k*g)*(exp(-k/m*t)-1)+g*t)+y_0;
          [~,idxMaxY] = max(y_pts);
          c(1) = -y_pts(idxMaxY)+ y_threshold;
          [~,idxYZero] = min(abs(y_pts(idxMaxY+1:end)));
          c(2) = -x_pts(idxMaxY+idxYZero) + x_threshold;
          f = m;
          ceq = [];
      end

      function stop = outfun(~,~,~)
          % Visualization
          stop = false;
          set(hopt,{'XData','YData'},{x_pts,y_pts});
      end

fprintf('v0=%f[m/s], theta=%f[deg], m=%f[kg]\n',x(1),rad2deg(x(2)),x(3));
end

軌道：

質量の変化：

%

[Shigenori Ishihara]

Tohru Kikawada様 有難うございます。早速提示頂いたサンプルを試してみました。 瑣末ですが、非線形不等号制約c(2)≦0を満たさずに終わっているのは　どうしてなのでしょうか？

[Tohru Kikawada]

c(2) の計算方法が正しくなかったので修正しました。

制約が守られていないなど終了したときの状態を知りたい場合には fmincon の戻り値の exitflag をご覧ください。