Problem writing fminsearch function

4 visualizaciones (últimos 30 días)
Wiktoria Schabek
Wiktoria Schabek el 3 de Jun. de 2022
Respondida: M Mirrashid el 5 de Jun. de 2022
I need help with code. I have identified the mathematical model as in the code below (J1,J2,J3). Each Ji includes a Pi in the range 2000-4000, but all J need to work together so all in all they work in P range 6000-12000. I needs to find the minimum value of the sum of all J. Below is my code that I try to write, but its give me error:
Error in optymalizacja2 (line 33) fun = @(P1,P2,P3) J (J1(P1)+J2(P2)+J3(P3));
Error in fminsearch (line 201) fv(:,1) = funfcn(x,varargin{:});
P1 = 2000 : 1 : 4000;
P2 = 2000 : 1 : 4000;
P3 = 2000 : 1 : 4000;
J1 = @(P1)(7.918*10.^(-20) * P1.^6) - (1.89*10.^(-15) * P1.^5) + (1.181*10.^(-11) * P1.^4) - (8.86*10.^(-8) * P1.^3) + (2.34*10.^(-4) * P1.^2) + (0.318 * P1) + 183.47 ;
J2 = @(P2)(5.41*10.^(-10) * P2.^3) - (3.51*10.^(-6) * P2.^2) + (0.071 * P2) + 3.69 ;
J3 = @(P3)(6.07*10.^(-9) * P3.^3) - (3.22*10.^(-5) * P3.^2) + (0.056 * P3)- 25.42 ;
fun = @(P1,P2,P3) J (J1(P1)+J2(P2)+J3(P3));
[x,fval] = fminsearch(fun, Rmin);
where:
R = @(P1,P2,P3)R1(P1)+R2(P2)+R3(P3);
% R is a function P - losses
Rmin = R(2000,2000,2000);
I think about use fmincon, but anyway I can't write it correctly.
  3 comentarios
Torsten
Torsten el 3 de Jun. de 2022
Editada: Torsten el 3 de Jun. de 2022
If you want to constrain P1,P2 and P3 as given, use fminsearchbnd:
fun = @(P) J1(P(1)) + J2(P(2)) + J3(P(3))
(I don't know what your J means in fun = @(P1,P2,P3) J (J1(P1)+J2(P2)+J3(P3));)
How R comes into play, I also could not understand.
Wiktoria Schabek
Wiktoria Schabek el 3 de Jun. de 2022
J in fun is accually my mistake, it should't be there R is function of output power which is dependent on input power (P1, P2 and P3) and function Li(Pi) (losses) Formula is: R1=P1-L1(P1) etc. Thanks for advise, I will try fminsearchbnd

Iniciar sesión para comentar.

Respuesta aceptada

Torsten
Torsten el 3 de Jun. de 2022
Editada: Torsten el 3 de Jun. de 2022
Here is a solution without optimizers:
% Define lower and upper bounds for P1, P2 and P3
lb1 = 2000;
ub1 = 4000;
lb2 = 2000;
ub2 = 4000;
lb3 = 2000;
ub3 = 4000;
% Define J1, J2 and J3 as polynomial functions
J1 = [7.918e-20 -1.89e-15 1.181e-11 -8.86e-8 2.34e-4 0.318 183.47];
J2 = [5.41e-10 -3.51e-6 0.071 3.69];
J3 = [6.07e-9 -3.22e-5 0.056 -25.42];
% Define dJ1/dP1, dJ2/dP2 and dJ3/dP3 as polynomial functions
dJ1 = [7.918e-20*6 -1.89e-15*5 1.181e-11*4 -8.86e-8*3 2.34e-4*2 0.318];
dJ2 = [5.41e-10*3 -3.51e-6*2 0.071];
dJ3 = [6.07e-9*3 -3.22e-5*2 0.056];
% Define d^2J1/d^2P1, d^2J2/d^2P2 and d^2J3/d^2P3 as polynomial functions
ddJ1 = [7.918e-20*6*5 -1.89e-15*5*4 1.181e-11*4*3 -8.86e-8*3*2 2.34e-4*2];
ddJ2 = [5.41e-10*3*2 -3.51e-6*2];
ddJ3 = [6.07e-9*3*2 -3.22e-5*2];
% Calulate roots of the derivatives
rJ1 = roots(dJ1)
rJ2 = roots(dJ2)
rJ3 = roots(dJ3)
% Select roots that could be minimum in the respective lb-ub intervals
idx1 = abs(imag(rJ1)) < eps & real(rJ1) >= lb1 & real(rJ1) <= ub1 & real(polyval(ddJ1,rJ1))>=0 ;
idx2 = abs(imag(rJ2)) < eps & real(rJ2) >= lb2 & real(rJ2) <= ub2 & real(polyval(ddJ2,rJ2))>=0 ;
idx3 = abs(imag(rJ3)) < eps & real(rJ3) >= lb3 & real(rJ3) <= ub3 & real(polyval(ddJ3,rJ3))>=0 ;
rJ1 = rJ1(idx1);
rJ2 = rJ2(idx2);
rJ3 = rJ3(idx3);
% Evaluate J1, J2 and J3 in the boundary points and in the interior points
% that could be minimum
x1 = [lb1,rJ1,ub1];
v1 = polyval(J1,x1);
x2 = [lb2,rJ2,ub2];
v2 = polyval(J2,x2);
x3 = [lb3,rJ3,ub3];
v3 = polyval(J3,x3);
[min1,idx1] = min(v1);
[min2,idx2] = min(v2);
[min3,idx3] = min(v3);
%Output P1_min, J1(P1_min), P2_min, J2(P2_min), P3_min and J3(P3_min)
x1(idx1)
min1
x2(idx2)
min2
x3(idx3)
min3
%Plot J1, J2, J3 over the lb-ub-intervals for comparison
x1 = lb1:1:ub1
plot(x1,polyval(J1,x1))
hold on
x2 = lb2:1:ub2
plot(x2,polyval(J2,x2))
hold on
x3 = lb3:1:ub3
plot(x3,polyval(J3,x3))
hold off

Más respuestas (1)

M Mirrashid
M Mirrashid el 5 de Jun. de 2022

Categorías

Más información sobre Econometrics Toolbox en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by