Series Summation With Function Handles

3 次查看（过去 30 天）
Andrew Poissant
Andrew Poissant 2017-4-27
评论： Andrew Newell 2017-4-27
I am performing a Monte Caro simulation with non-linear optimization inside using fmincon. Fmincon requires the use of function handles. However, I want to perform a series summation, as attempted when defining Elife with symsum. Since symsum only takes symbolic variables, how would I perform this summation using function handles? To simplify the code, everything before the for loop is just setting up the distributions and inside the for loop is where I define a couple functions and where I set up and perform the optimization.
clear all
m_tlife = 33;
std_tlife = 11;
dist_tlife = makedist('Normal', m_tlife, std_tlife);
a_drate = 2;
b_drate = 0.006;
dist_drate = makedist('Gamma', a_drate, b_drate);
m_Ee = 4.56;
std_Ee = 0.27;
dist_Ee = makedist('Normal', m_Ee, std_Ee);
m = 34;
Asp = 40.1; % m^2
em = 0.162;
Np = 10;
for i = 1:Np
tlife = random(dist_tlife);
sigma = random(dist_drate);
Ee = random(dist_Ee);
Ep = Asp*em*Ee*365/m;
Elife = @(x)symsum(Ep*(m - x(2))*(1 - sigma)^x(1) + Ep*x(2)*(1 - sigma)^(tlife - x(1)), x(1), 1, tlife);
fun = @(x)Elife(x)*x(2)/(1 + x(1))
x0 = [0,0];
lb = [0,0];
up = [tlife m];
A = [];
b = [];
Aeq = [];
beq = [];
x(i,:) = fmincon(fun, x0, A, b, Aeq, beq, lb, up);
L(i) = fun(x)*10^2;
end

请先登录，再回答此问题。

回答（1 个）

Andrew Newell
Andrew Newell 2017-4-27
In principle, you can use matlabFunction to get around this difficulty. However, I noticed a couple of issues. First, the summation is over a variable tlife that is generally not an integer, which is an invalid choice for the last parameter in symsum. So you may have to round it off.
Second, to define Elife, you need symbolic variables to stand in for x(2) and x(1). I tried the following (for reproducibility, I hard coded the parameter values in):
tlife = 28.2305;
sigma = 0.0712;
Ee = 5.3077;
Ep = 370.1550;
syms x1 x2
Elife = symsum(Ep*(m - x2)*(1 - sigma)^x1 + Ep*x2*(1 - sigma)^(tlife - x1), x1, 1, round(tlife));
This gave the result
(375669370696008546723677936360824844289937144399794339030368775708855024128885589600699801*1161^(461/2000)*1250^(1539/2000)*x2)/103397576569128459358926086508745356695726513862609863281250000000000000000000000000000000 - (436152139378065922746190084114917644220617024648161227614258148597980683013636169526412468961*x2)/103397576569128459358926086508745356695726513862609863281250000000000000000000000000000000 + 7414586369427120686685231429953599951750489419018740869442388526165671611231814881949011972337/51698788284564229679463043254372678347863256931304931640625000000000000000000000000000000
which is surprising because x1 has disappeared. Thus, in applying matlabFunction, I get a function only of the second component:
Elife = matlabFunction(symsum(Ep*(m - x2)*(1 - sigma)^x1 + Ep*x2*(1 - sigma)^(tlife - x1), x1, 1, round(tlife)));
disp(Elife)
@(x2)x2.*(-4.218204660594419e3)+5.086790208514118.*2.415862736086788e2.*x2.*3.633251214982273+1.434189584602102e5
Given this result, the only way I could find to define the function fun that didn't give an error called the second component explicitly:
fun = @(x) Elife(x(2))*x(2)/(1+x(1));
This does get minimized without error, hitting one of the constraints.
In summary, I ran the whole loop successfully using the following block of code instead of your definition of fun:
syms x1 x2
Elife = matlabFunction(symsum(Ep*(m - x2)*(1 - sigma)^x1 + Ep*x2*(1 - sigma)^(tlife - x1), x1, 1, round(tlife)));
fun = @(x) Elife(x(2))*x(2)/(1+x(1));
But if I were you, I'd try to understand why x1 is disappearing from Elife.
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2017-4-27
I agree, after the symsum the index variable would be expected to disappear.
Andrew Newell
Andrew Newell 2017-4-27
Oh, right! Not seeing the forest for the trees.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by