is it possible to determine the values ​​of the input parameters ca, cb, cc and cd for which the cc value at the output would be maximal?

is it possible to determine the values ​​of the input parameters ca, cb, cc and cd for which the cc value at the output would be maximal?

    dca/dt=-theta(1) * ca-theta(2).cA;
    dcb/dt= theta(1).*cA+theta(4).*cc*theta(3).*cb-theta(5).*cd
    dcc/dt= theta(2).*cA+theta(3).*cb*theta(4).*cCtheta(6).*cd
    dd/dt= theta(5).*cb*theta(6).*cd

Values for all thetas are known.

Thanks in advance.

[Walter Roberson]

I think we need a correction near cC in dcc/dt

Is the last line dcd/dt instead of dd/dt ?

is CA a constant or a typo?

[Walter Roberson]

The dcc/dt tells us that cc is a function of t. It might not be possible to find functions ca, cb, cc, cd such that cc(t) is maximal for all t, compared to all of the other possible ca, cb, cc, cd functions. It would make more sense to talk about maximal at a particular t.

[Walter Roberson]

Assuming that the third line is cC - theta(6) and that cA and cC are constant and that the last line is dcd/dt, then there are three forms of solution for the differential equation. Two of them have to essentially be defined point by point. The third one has cb(t) = 0 and leads to a cc(t) that is linear in t, with cd(t) being constant, and ca(t) having an exponential term.

[Bosnian Kingdom]

I think we need a correction near cC in dcc/dt

I made mistake with dcc/dt, the corect is dc/dt.

Is the last line dcd/dt instead of dd/dt ?

No, it is dd/dt.

is CA a constant or a typo?

No, CA is ca and is changed depending on t.

So the correct set of equations are:

dca/dt=-theta1 * ca-theta2*ca;
dcb/dt= theta1*ca+theta*cc*theta3*cb-theta5*cd
dcc/dt= theta2*ca+theta3*cb*theta4*cc*theta6*cd
dd/dt= theta5*cb*theta6*cd

Thetas are known (fixed values)

t is also known (t at the beggining and t finally).

Is it possible to find inlet values for ca, cb, cc ans cd for the maximal value of finally cc ?

[Torsten]

If it's "dd/dt = theta5*cb*theta6*cd", where do you get "cd" from ?

And concerning your question about maximal value of final cc:

Choose initial values ca0 = cb0 = cd0 = 0 and cc0 = +Infinity.

[Bosnian Kingdom]

You are right:

dd is actually dcd, so the equatons are:

dca/dt=-theta1 * ca-theta2*ca;
dcb/dt= theta1*ca+theta*cc*theta3*cb-theta5*cd
dcc/dt= theta2*ca+theta3*cb*theta4*cc*theta6*cd
dd/dt= theta5*cb*theta6*cd

Is your answer still:

Choose initial values ca0 = cb0 = cd0 = 0 and cc0 = +Infinity.

[Torsten]

dcd/dt= theta5*cb*theta6*cd

And yes, the answer still remains

Choose initial values ca0 = cb0 = cd0 = 0 and cc0 = +Infinity

unless there are no restrictions on the initial concentrations of the species.

[Bosnian Kingdom]

I tried, and whatever initial values out i get the same results. Is there some objection function which i can write to maximize Cc? Thanks for answer above.

[Torsten]

Without seeing your code, nobody will be able to help.

[Bosnian Kingdom]

You are right.

First I made function:

function diffeqs = ode(t,var)

ca=var(1);
cb=var(2);
cc=var(3);
cd=var(4);
theta1 =  0.89197;
theta2 =  0.26069;
theta3 =  0.24870;
theta4 =  0.48183;
theta5 =  0.60977;
theta6 =  0.01301;

diffeqs(1,1)=-theta1 * ca-theta2*ca;
diffeqs(2,1)= theta1*ca+theta4*cc*theta3*cb-theta5*cd;
diffeqs(3,1)=theta2*ca+theta3*cb*theta4*cc*theta6*cd;
diffeqs(4,1)=theta5*cb*theta6*cd;

end

Then I made script:

range=[0:6];
ICs=[0,0,NaN,0];

[tsol,varsol]=ode45(@ode,range,ICs);

And when I run, whatever initial values i set, the result are the same.

[Walter Roberson]

That NaN is sure to cause problems.

There is a solution to the system of

ca(t) = 0, cb(t) = 0, cc(t) = cc0, cd(t) = 0

where cc0 is the initial condition for cc. WIth you having provided 0 as the initial conditions for ca, cb, cd, then you would possibly be hitting this situation.

There appears to be another solution when cb(0) is non-zero:

ca(t) = 0, cb(t) = cb0 * exp(theta3*theta4*cc0*t), cc(t) = cc0, cd(t) = 0

When cb(0) = 0, then there is a solution with a non-zero ca with a simple form, and with cd(t) = 0, and with a complicated cb expression, and with

cc(t) = (-exp(-(theta1+theta2)*t) * ca0 * theta2 + (ca0+cc0) * theta2 + cc0 * theta1) / (theta1 + theta2)

That system is linear in cc0, so you would maximize it when cc(0) is infinity. (The system has cc0*(theta1+theta2) in the numerator, and has (theta1+theta2) in the denominator, so if theta1+theta2 is negative, the two negatives would cancel out as cc0 reached infinity, and so the maxima is at cc(0) -> infinity rather than at cc(0) -> infinity * sign(theta1+theta2) )

The complicated expression for cb does not have a closed form (except perhaps at very specific reationships between the angles), but that does not matter since you just want to know for cb, which does have a simple relationship in this case.

When cc(0) = 0, then there is an additional fairly messy solution that cannot be expressed in closed form, and which is difficult to to maximize. Oddly, even though the expression involves cb0 in non-trivial ways, the derivative is 0 with respect to cb0, suggesting that perhaps the complicated expression could be simplified... just maybe to a fairly simple expression ??

[Bosnian Kingdom]

Thanks for you answer.

Can you explain me, where i will write:

ca(t) = 0, cb(t) = 0, cc(t) = cc0, cd(t) = 0


ca(t) = 0, cb(t) = cb0 * exp(theta3*theta4*cc0*t), cc(t) = cc0, cd(t) = 0


or 

cc(t) = (-exp(-(theta1+theta2)*t) * ca0 * theta2 + (ca0+cc0) * theta2 + cc0 * theta1) / (theta1 + theta2)

Will i write this in function or script?

[Walter Roberson]

function [ca, cb, cc, cd] = generate1(theta, ca0, cb0, cc0, cd0)
   ca = @(t) zeros(size(t));
   cb = @(t) zeros(size(t));
   cc = @(t) cc0 .* ones(size(t));
   cd = @(t) zeros(size(t));
end

function [ca, cb, cc, cd] = generate2(theta, ca0, cb0, cc0, cd0)
  ca = @(t) zeros(size(t));
  cb = @(t) cb0 .* exp(theta(3) .* theta(4) .* cc0 .* t);
  cc = @(t) cc0 .* ones(size(t));
  cd = @(t) zeros(size(t));
end

function [ca, cb, cc, cd] = generate3(theta, ca0, cb0, cc0, cd0)
  ca = @(t) error('sorry, not available at the moment, I closed that page');
  cc = @(t) (-exp(-(theta(1) + theta(2)).*t) .* ca0 .* theta(2) + (ca0+cc0) .* theta2 + cc0 .* theta(1)) ./ (theta(1) + theta(2));
  cd = @(t) zeros(size(t));
  error('sorry, cb is pretty messy involving an integral, and I closed that page');
end

[Bosnian Kingdom]

Thank you for the answer:

I do this:

function diffeqs = ode(t,var)

ca=var(1);
cb=var(2);
cc=var(3);
cd=var(4);
function [ca, cb, cc, cd] = generate1(theta, ca0, cb0, cc0, cd0);
   ca = @(t) zeros(size(t));
   cb = @(t) zeros(size(t));
   cc = @(t) cc0 .* ones(size(t));
   cd = @(t) zeros(size(t));
end
theta1 =  0.89197;
theta2 =  0.26069;
theta3 =  0.24870;
theta4 =  0.48183;
theta5 =  0.60977;
theta6 =  0.01301;

diffeqs(1,1)=-theta1 * ca-theta2*ca;
diffeqs(2,1)= theta1*ca+theta4*cc*theta3*cb-theta5*cd;
diffeqs(3,1)=theta2*ca+theta3*cb*theta4*cc*theta6*cd;
diffeqs(4,1)=theta5*cb*theta6*cd;

end

and

range=[0:6];
ICs=[1,0,0,0];

[tsol,varsol]=ode45(@ode,range,ICs);

And I got this results:

1 0 0 0

0.315743172332825 0.533782793523100 0.154754144678011 0

0.0997404581818380 0.715101440379562 0.203606145746861 0

0.0315030140235672 0.787188633130357 0.219038987450069 0

0.00994818871230687 0.825351252258936 0.223913909292062 0

0.00314254428025749 0.853241325366288 0.225453099900734 0

0.000992588060197626 0.878328022861088 0.225939342233258 0

What is the problem?

[Walter Roberson]

You do not define ca0 cb0 cc0 cd0

[Bosnian Kingdom]

I define ca0, cb0, cc0 and cd0. And i still getting same results:

function diffeqs = ode(t,var)
ca=var(1);
cb=var(2);
cc=var(3);
cd=var(4);
function [ca, cb, cc, cd] = generate1(theta, ca0, cb0, cc0, cd0)
    ca0=0;
    cb0=0;
    cc0=0;
    cd0=0;
   ca = @(t) zeros(size(t));
   cb = @(t) zeros(size(t));
   cc = @(t) cc0 .* ones(size(t));
   cd = @(t) zeros(size(t));
end
theta1 =  0.89197;
theta2 =  0.26069;
theta3 =  0.24870;
theta4 =  0.48183;
theta5 =  0.60977;
theta6 =  0.01301;
diffeqs(1,1)=-theta1 * ca-theta2*ca;
diffeqs(2,1)= theta1*ca+theta4*cc*theta3*cb-theta5*cd;
diffeqs(3,1)=theta2*ca+theta3*cb*theta4*cc*theta6*cd;
diffeqs(4,1)=theta5*cb*theta6*cd;
end

[Walter Roberson]

What results are you expecting? With initial conditions 0 your equations only have the all 0 solution.

[Bosnian Kingdom]

I don't know initial values, i need to determine initial values for maximal output concentration cc. Is it posible?

[Walter Roberson]

I went through the scenarios in

https://www.mathworks.com/matlabcentral/answers/458387-is-it-possible-to-determine-the-values-of-the-input-parameters-ca-cb-cc-and-cd-for-which-the-cc#comment_701095

If you start with ca 0, cb 0, cd 0, then the output for cc will be the same as your input for cc. Therefore you can get any output you want for cc by using that same input for cc. This is the same as Torsten told you, that to maximize cc then start with cc being +infinity.

This is the "generate1" scenario. You start with all 0 and you get the same cc output as input. Nothing ever changes.

The "generate2" scenario starts with a non-zero cb, and has a cb that changes over time. However, the cc stays the same as input, so to maximize cc, start with cc being infinity.

The "generate3" scenario does have a cc that changes. But look at it: it has terms that are independent of time, and it has a term which is -ca0*theta2*exp(-(theta1+theta2)*t) / (theta1+theta2) .

If theta1+theta2 is positive, then -(theta1+theta2)*t gets more negative as time increases, and exp() of a more negative term decreases towards 0, so -ca0*theta2 times the term will get closer and closer to 0 (even if ca0 or theta2 are negative.) Thus in that situation, yes, cc would increase, but it would increase towards a constant upper bound determined by the terms that are independent of time.

If theta1+theta2 is negative, then -(theta1+theta2)*t gets more positve as time increases, and exp() of a more positive term increases rapidly, and the negative of that grows rapidly smaller. But with the division by (theta1+theta2) and (theta1+theta2) being negative, that would give you a term that would grow rapidly larger in time. But theta1+theta2 would have to be negative, and your input theta are all positive.

If you have to maximize cc, and you want one of the situations where cc is not constant, then you are down to two possibilities:

1. The generate3 scenario where cc starts as small as it can ever get and grows towards a constant upper bound as the exponential contribution of it decreases until you can ignore it; or
2. The formulation with ca(0) non-zero, cb(0) non-zero, cc(0) non-zero, cd(0) irrelevant, in which case the situation is too complicated to analyze analytically without a lot of work. ca(t), cb(t), cc(t) would all have to be expressed in terms of functions of integrals.

[Bosnian Kingdom]

How can I accept you answer? You helped me. I used generate3. Thank you.