How do I optimize solutions of ODE?

I'm trying to follow the process below. I would appreciate if you could let me know how to do this.

① solve the ODE which includes an optimization variable.

② optimize the objective function, which is the solution of ODE at a particular point.

I also post the code I wrote. I think the ODE needs to be converted in some way, but I don't know how.

c = optimvar("c","LowerBound",0,"UpperBound",100); 
Fc = 10; 


M = @(t,Y)[Y(2);sign(Y(2)).*(-1.0./3.0e+1).*Fc-Y(1).*1.048982544e+3-Y(2)./3.0e+1.*c]; %ODE


interval=[0:0.001:0.5];

yinit=[0.008337,0];

[t,y] = ode45(M,interval,yinit); 

Error using superiorfloat (line 13)
Inputs must be floats, namely single or double.

Error in odearguments (line 114)
dataType = superiorfloat(t0,y0,f0);

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

prob = optimproblem;

prob.Objective = 0.00592708-y(1,194); %objective function

c0.c = 0; 

sol = solve(prob,c0)

Answer 1

[Sam Chak]

Hi @Amito Usami-san

Considering your objective of varying the force input 'Fc', it becomes necessary to create a Gain Schedule for the damping coefficient 'c' as a function of the force. Through a simple trial-and-error approach, it appears that a linear function for the damping coefficient c can yield relatively satisfactory results.

Please take a look at the Gain Schedule provided below:

format long

%% Gain Schedule for damping coefficient c
c       = @(Fc) 105.030559597794 - 5.54966340482472*Fc;
force   = linspace(5, 15, 201);
plot(force, c(force)), grid on
xlabel('Force'), ylabel('c'), title('Gain Schedule for damping coefficient c')

%% 2nd-order ODE
Fc      = 8;            % Force input
M       = @(t,Y) [Y(2);
                  sign(Y(2))*(- 1.0/3.0e+1)*Fc - Y(1)*1.048982544e+3 - Y(2)/3.0e+1*c(Fc)];

%% call ode45 solver
interval= [0:0.001:0.5];
yinit   = [0.008337, 0];
[t, y]  = ode45(M, interval, yinit); 

%% plot results
out     = y(:,1);
[pks, locs] = findpeaks(out)

pks = 2x1

0.005926983356573 0.003944587078435

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

locs = 2x1

195 390

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

idx     = locs(1);
plot(t, out), hold on
plot(t(idx), out(idx), 'o', 'markersize', 12), grid on
xlabel('t'), ylabel('y(t)'), title('Response and Desired Peak')
legend('Response', 'Desired Peak')

[Amito Usami]

Sorry if I didn't explain it well enough. What I'm trying to do is to write a code to find the damping-coefficient c such that the value of y(1,194) remains the same as 0.00592708 if the value of Fc changes. For example, when Fc = 9, c would be about 54.9, and when Fc = 10, c would be about 49.4.

[Sam Chak]

No worries. What does y(1,194) mean? Is it a specific location in the y vector?

[Amito Usami]

Sorry again, it is not y(1,194) but y(194,1). y(194,1) means the second largest amplitude when I solve this ODE normally below.

Fc = 10;
c = 49.4;


M = @(t,Y)[Y(2);sign(Y(2)).*(-1.0./3.0e+1).*Fc-Y(1).*1.048982544e+3-Y(2)./3.0e+1.*c];
interval=[0:0.001:0.5];
yinit=[0.008337,0];
[t,y]=ode45(M,interval,yinit);

plot(t,y(:,1))

title("displacement response")
xlabel("time(s)")
ylabel("displacement(m)")
ylim([-0.01 0.01])
grid on
ax=gca;
ax.XAxisLocation='origin';
ax.FontSize=10;
ax.TitleFontSizeMultiplier=1.5;

y(194,1)

ans = 0.0059

[Sam Chak]

You want the second highest peak to be exactly 0.00592708?

[Amito Usami]

yes.

[Sam Chak]

Hi @Amito Usami

I understand how you varied the value for 'c' to achieve the desired peak. Based on that, I have made updates to the answer by implementing a Gain Schedule for the damping coefficient 'c' as a function of the force input 'Fc'.

Please let me know if the solution provided by the Gain Schedule is satisfactory to you.

[Amito Usami]

Thank you. Can I ask you how to get this equation?

c = @(Fc) 105.030559597794 - 5.54966340482472*Fc;

[Sam Chak]

Hi @Amito Usami

Similar to what you did for (Fc = 9, c = 54.9) and (Fc = 10, c = 49.4), I conducted manual tuning of the parameter c across a wider range of values for Fc. Through this process, I discovered a linear trend. Consequently, it became straightforward to determine the equation of the straight line, .

That's why it took me a few hours to reply to you. If you find the solution helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Your support is greatly appreciated!

Here is just a demo using your original two-point data:

Fc      = [9, 10];
c       = [54.9, 49.4];
csch    = fit(Fc', c', 'poly1') % gain schedule for c

csch = 
     Linear model Poly1:
     csch(x) = p1*x + p2
     Coefficients:
       p1 =        -5.5
       p2 =       104.4

[Amito Usami]

I've got it. Thank you so much.

[Sam Chak]

You are welcome, @Amito Usami. It's important to note that this approach is known as Gain-scheduling. It is quite similar to creating a lookup table with an adequate number of data points and subsequently interpolating the data using the curve-fitting method.