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

Global Optimization ToolboxoptimizationOptimization Toolboxsimulation日本語

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

Best Answer

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
軌道：
質量の変化：
%
Related Question