Effacer les filtres
Effacer les filtres

Problem writing fminsearch function

4 vues (au cours des 30 derniers jours)
Wiktoria Schabek
Wiktoria Schabek le 3 Juin 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 commentaires
Torsten
Torsten le 3 Juin 2022
Modifié(e) : Torsten le 3 Juin 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 le 3 Juin 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

Connectez-vous pour commenter.

Réponse acceptée

Torsten
Torsten le 3 Juin 2022
Modifié(e) : Torsten le 3 Juin 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

Plus de réponses (1)

M Mirrashid
M Mirrashid le 5 Juin 2022

Catégories

En savoir plus sur Econometrics Toolbox dans Help Center et File Exchange

Produits


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by