0 个投票

Hey everyone,
I want to make long-term load forecasting using GA. So the first step is to come up with a model, in one of the papers the objective function is a polynomial of tenth order:
obj= c10*x.^10 + c9*x.^9 + c8*x.^8 + c7*x.^7 + c6*x.^6 + c5*x.^5 + c4*x.^4 + c3*x.^3 + c2*x.^2 + c1*x.^1 + c0*x.^0;
In order to make the obj function ready for the GA I need to estimate the coefficients.
The rest of my code is as follows:
>> f = @(c,x) 1 + c(1)*x.^1 + c(2)*x.^2 + c(3)*x.^3 + c(4)*x.^4 + c(5)*x.^5 + c(6)*x.^6 + c(7)*x.^7 + c(8)*x.^8 + c(9)*x.^9 + c(10)*x.^10;
>> cfit = nlinfit(xdata,ydata,f,c)
all the data that I have are the years from 1982 till 1991 and the corresponding demand in each year.
I didn't understand nlinfit quite well,, what I am supposed to put in place of xdata, ydata and c.
Any help will be appreciated.

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

 采纳的回答

Image Analyst
Image Analyst 2011-12-19

2 个投票

Why do you call that non-linear regression? It's just a regular polynomial and it's linear in the coefficients, c. You don't have c(6)^2 or log(c(5)) or anything non-linear like that. It's just c(#) to the first power multiplied by the x to some power. Because your x are non-linear does not make it non-linear regression. All your c's are linear so it's linear regression. So you can simply use polyfit() and simplify your life.

9 个评论
显示 7更早的评论 隐藏 7更早的评论

Yasmin Tamimi
Yasmin Tamimi 2011-12-19
OK thnx a lot,,
so shall I write y as the years in an ascending or descending order or it won't differ?!!
bym
bym 2011-12-19
one wonders why somebody needs a 10th order polynomial <8-0
Yasmin Tamimi
Yasmin Tamimi 2011-12-19
I had the same inquiry, but unfortunately the paper I am reading from didn't include its derivation!! I tried to contact the author but no reply!!
Image Analyst
Image Analyst 2011-12-19
Unfortunately MATLAB does it in reverse of the way you'd think and in reverse of the way I usually see it done. c(1) is the coefficient of the highest order x term, x^n, while c(n) is the coefficient of x (not x^n) and c(n+1) is the additive constant term.
Image Analyst
Image Analyst 2011-12-19
I'd say that you can NOT use a 10th order term for forecasting. You'll have a heck of a time getting valid values just within the training zone, let alone forecasting (extrapolating)!!! A 10th order polynomial will oscillate wildly (even within the training zone) and rapidly zoom off to + or -infinity once you start extrapolating, basically giving garbage there.
Walter Roberson
Walter Roberson 2011-12-19
A 10th order polynomial will have 10 maxima, 10 minima and 10 zeros, which _I_ would not consider to be "oscillating wildly". But it will indeed rapidly go off the + or - infinity.
Image Analyst
Image Analyst 2011-12-20
In examples I've seen, the points "in between" two "training" points tend to move farther away from a straight line in between the two training points as the order of the equation grows larger. That's what I meant. And those in between points are a lot more sensitive to slight changes in the training points locations.
Walter Roberson
Walter Roberson 2011-12-20
Hmmm, that's probably provable, too -- though it could plausibly be the case that if the zeros were carefully positioned at decreasing intervals that at least one measure of the swing might decrease.... Yes, indeed, I have just constructed a sequence whose zeros do not change, but whose maximums swing less and less as the length of the sequence increases. Certainly, though, in my first trial series the maximums increased distinctly as the length of the sequence increased.
Can I know are there any other methods I can use also to compare the performances among methods. I used nlinfit and lsqcurvefit, I looked up and found fitnlm and lsqnonlin are same as the about methods. And I have looked into different methods such as ridge , robust , polyfit but none of them fit the case that lsqcurvefit is considering : as in lsqcurvefit(fun,x0,xdata,ydata) *nonlinear case Please help me =( , I have been looking at it for a while

请先登录，再进行评论。

更多回答（5 个）

Richard Willey
Richard Willey 2011-12-19

1 个投票

I'd strongly suggestion that you watch a webinar titled "Electricity Load and Price Forecasting with MATLAB". The webinar is available at: http://www.mathworks.com/company/events/webinars/wbnr51423.html
All of the code and the data sets are available on MATLAB Central.
This webinar shows two different ways to model the demand for electric power. The first is based on a neural network. The second uses bagged decision trees. The code also includes safeguards to protect against overfitting.
I'm also going to point you at a blog posting that I wrote on data driven fitting. If you are primarily worried about interpolation you might find this a useful alternative to high order polynomials
http://blogs.mathworks.com/loren/2011/01/13/data-driven-fitting/

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Yasmin Tamimi
Yasmin Tamimi 2011-12-19
Thnx a lot Richard for the webinar, but I have to use GA instead of NN and long-term forecasting instead of short.

请先登录，再进行评论。

Greg Heath
Greg Heath 2011-12-19

0 个投票

You definitely do not want a high order polynomial for prediction.
Check out Richard's references.
Greg
Yasmin Tamimi
Yasmin Tamimi 2011-12-19

0 个投票

Actually the prediction went wrong, it seems like I am having an error during the running of GA!! so I minimized the order till 2 and I still have the same error!! here's my code:
FIRST M-FILE:
format long e
f = @(c,x) c(1)*x^2 + c(2)*x^1 + c(3)*x^0;
% I have 20 data points for both the years and the load but I should use the first 10 to calculate the coefficients and the other 10 should be predicted using ga:
years = [1982 1983 1984 1985 1986 1987 1988 1989 1990 1991];
load = [1702 2344 2097 2313 2588 2885 4341 4779 5251 5721];
estimated_coefficients = polyfit(years,load,2);
% c1 = 4.165909091390513e+001;
% c2 = -1.650490772918446e+005;
% c3 = 1.634786853371606e+008;
SECOND M-FILE:
%% The objective function function y = load_forecast(x)
c1 = 4.165909091390513e+001;
c2 = -1.650490772918446e+005;
c3 = 1.634786853371606e+008;
y = c1*x(1)^2 + c2*x(2)^1 + c3;
THIRD M-FILE:
FitnessFcn=@load_forecast;
GenomeLength = 2; % Number of variables in the fitness function
LB = zeros(1,2); % Lower bound
UB = ones(1,2); % Upper bound
Bound = [LB;UB];
% options structure
options = gaoptimset('Vectorized','on','PopulationType','bitstring','CreationFcn',@int_pop,'MutationFcn',{@mutationuniform,0.04},... 'CrossoverFcn',{@crossoverscattered,0.8}, 'PopInitRange' ,Bound, 'Display','Iter','StallGenL',100,'Generations',150, ... 'PopulationSize',50);
[X,FVAL] =ga(@load_forecast,2,[],[],[],[],LB,UB,[],options);
AND THE ERROR THAT I GET IS:
??? Reference to non-existent field 'Verbosity'. Error in ==> gacommon at 79 [Iterate.x,Aineq,bineq,Aeq,beq,lb,ub,msg,exitFlag] = ... Error in ==> ga at 269 [x,fval,exitFlag,output,population,scores,FitnessFcn,nvars,Aineq,bineq,Aeq,beq,lb,ub, ... Error in ==> ga_load_forecast at 27 [X,FVAL] =ga(@load_forecast,2,[],[],[],[],LB,UB,[],options);
FINAL QUESTION: the data that I have shouldn't it be incorporated
within the fitness function or GA in any way or another??
really any help is appreciated..

12 个评论
显示 10更早的评论 隐藏 10更早的评论

bym
bym 2011-12-20
I don't know anything much about GA but you can use
y = polyval(estimate_coefficients,x)
for your 2nd m-file. Also, shouldn't the objective function be some measure of the error (residual); you are just returning the function value ?!?
Yasmin Tamimi
Yasmin Tamimi 2011-12-20
You mean after using y = polyval(estimate_coefficients,x)
this will be my fitness function??
I want to use ga so that the I can get the load demand which is so close to the real value.
Yes, the error was also one of my concerns, till now i couldn't figure it out bcz I am getting errors!!
Image Analyst
Image Analyst 2011-12-20
I wouldn't call it a "fitness" function. Those are your "fitted" or, more properly, "estimated" values.
Yasmin Tamimi
Yasmin Tamimi 2011-12-20
When I run polyfit command then polyval for the estimated coefficients,, the estimated load values are far away different from the actual ones and I am wondering why?!
Walter Roberson
Walter Roberson 2011-12-20
"??? Reference to non-existent field 'Verbosity'" -- so what happens if you add a Verbosity option in your gaoptimset ?
Greg Heath
Greg Heath 2011-12-23
Is
y = c1*x(1)^2 + c2*x(2)^1 + c3;
a misprint?
Greg
Walter Roberson
Walter Roberson 2011-12-23
It looks consistent with the other code to me?
Yasmin Tamimi
Yasmin Tamimi 2011-12-23
@ walter: actually i don't even know what verbosity means!
Yasmin Tamimi
Yasmin Tamimi 2011-12-23
@ greg: no it is not a misprint, i'm having a second order polynomial and i want to calculate its coefficients by using ga, here's the final code that i wrote:
GenomeLength = 3; % Number of variables in the fitness function: 3 coefficient
% options structure
options = gaoptimset('MutationFcn',{@mutationuniform,0.04},...
'CrossoverFcn',{@crossoverintermediate,0.8},'Display','iter','StallGenL',40,'Generations',50, ...
'PopulationSize',100);
x =ga(@load_forecast,3,[],[],[],[],[],[],[],options)
fid= fopen('Load Forecast using GA.txt','w');
fprintf(fid,'Years\t Actual Load\t Estimated Load\t Percentage Error\r\n');
fprintf(fid,'________________________________________________________________\r\n') ;
Years2 = [1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001];
Actual_Load2 = [1702 2344 2097 2313 2588 2885 4341 4779 5251 5721 6322 6992 7738 8570 9397 10338 11371 12507 13759 15134];
Estimated_Load2 = polyval(x,Years2);
Residual2 = Actual_Load2-Estimated_Load2;
Percentage_Error = (abs((Residual2)./ Actual_Load2)).*100;
for i = 1:length(Years2)
fprintf(fid,'%d\t %d\t\t %d\t %d\r\n',Years2(i),Actual_Load2(i),Estimated_Load2(i),Percentage_Error(i));
fprintf(fid, '\r\n');
end
plot(Years2,Actual_Load2,'bo-')
hold on
plot(Years2,Estimated_Load2,'rx-')
hold off
title('Actual Vs Estimated Load/Demand');
xlabel('Time (Years)');
ylabel('Load / Demand (MWatt)');
legend('Actual Load','Forecasted Load');
grid on
fclose(fid);
Yasmin Tamimi
Yasmin Tamimi 2011-12-23
yes although the code seems correct but it gave me answers far away from the actual!!
Yasmin Tamimi
Yasmin Tamimi 2011-12-23
here is the link to the paper i was trying to simulate part of their findings (2nd order poly. using their fitness function):
but i used the data from another paper bcz i didn't have access to reference [6] where they got the data from:
http://www.waset.org/journals/waset/v6/v6-32.pdf
i hope this will help you understand what i was doing.
Yasmin Tamimi
Yasmin Tamimi 2011-12-23
sorry i forgot to post my ff, here it is:
function y = load_forecast(x)
k = 0.0001; % k is a scaling constant
Actual_Load1 = [1702 2344 2097 2313 2588 2885 4341 4779 5251 5721];
T = [1982 1983 1984 1985 1986 1987 1988 1989 1990 1991]; % Years
t = 10; % number of years
% Here x is considered as the coefficient of the second order polynomial equal,
% and we want to find their optimal values such that the error(residual) is minimized
sum = 0;
for i = 1:t
sum_residual = abs(((x(1)*(T(i))^2) + (x(2)*T(i)) + x(3))- Actual_Load1(i));
sum = sum_residual;
end
y = 1 + (k * sum_residual);

请先登录，再进行评论。

Richard Willey
Richard Willey 2011-12-20

0 个投票

For what its worth, I just took a very quick look at the data set that you provided.
years = [1982 1983 1984 1985 1986 1987 1988 1989 1990 1991];
load = [1702 2344 2097 2313 2588 2885 4341 4779 5251 5721];
You can fit the years 1988 --> 1991 with an almost perfectly straight line. In a similar fashion, the years 1984 --> 1987 with another straight line. In both cases the R^2 is over .995.
I really don't understand that approach that you're taking... I feel like you're trying to force Genetic Algorithms into the solution space regardless of whether this is warranted.
Given that you're primarily interested in using GA, there's one last resource that I'd recommend looking at:
The "Global Optimization with MATLAB Products" provides a very good introduction to GA. You can watch the webinar at: http://www.mathworks.com/company/events/webinars/wbnr43346.html?seq=1
All of the code is available for download from MATLAB Central.

2 个评论
显示 无 隐藏 无

Yasmin Tamimi
Yasmin Tamimi 2011-12-20
Thnx a lot for the webinar. I am from the part of estimating the coefficients and building my model.
My only problem now is in writing the fitness function for the ga!!
Alex
Alex 2012-9-25
Yasmine, can you solve this? What do you do?

请先登录，再进行评论。

anahita
anahita 2025-3-24

0 个投票

hello, i want to write a script of nonelinear regression but i have a problem to write that because first idont know which one works for me so i wrote a code for nonelinear least squared regression but th R was 0.56. caan you help me in this clc; clear; close all;
data = readmatrix('analysis_1.xlsx', 'Sheet', 'Main');
S = data(:, 6);
mu = data(:, 2);
sigma = data(:, 3);
gamma = data(:, 4);
Y = data(:, 5);
ZR = data(:, 7);
%% Ryan-Joiner mothode
sorted_S = sort(S);
n = length(S);
z_scores = norminv(((1:n) - 0.5) / n, 0, 1);
S_mean = mean(S);
S_std = std(S);
e_i = (sorted_S - S_mean) / S_std;
Rp = sum(e_i .* z_scores) / sqrt(S_std^2 * (n - 1) * sum(z_scores.^2));
disp(['Ryan-Joiner Rp = ', num2str(Rp)]);
%% (Nonlinear Least Squares)
% S = a * mu^b * sigma^c * gamma^d * Y^e * ZR^f
modelFun = @(params, X) params(1) * (X(:,1) .^ params(2)) .* ...
(X(:,2) .^ params(3)) .* (X(:,3) .^ params(4)) .* ...
(X(:,4) .^ params(5)) .* (X(:,5) .^ params(6));
X = [mu, sigma, gamma, Y, ZR];
initial_guess = [1, 1, 1, 1, 1, 1];
params_opt = lsqcurvefit(modelFun, initial_guess, X, S);
a = params_opt(1);
b_coeff = params_opt(2:end);
disp(['a = ', num2str(a)]);
disp(['b = ', num2str(b_coeff(1))]);
disp(['c = ', num2str(b_coeff(2))]);
disp(['d = ', num2str(b_coeff(3))]);
disp(['e = ', num2str(b_coeff(4))]);
disp(['f = ', num2str(b_coeff(5))]);
%% R²
S_Estimate = modelFun(params_opt, X);
SS_res = sum((S - S_Estimate).^2);
SS_tot = sum((S - mean(S)).^2);
R2 = 1 - (SS_res / SS_tot);
disp(['R2 = ', num2str(R2)]);
results_table = table(S, S_Estimate, 'VariableNames', {'Actual_Storage', 'Estimated_Storage'});
coeff_table = array2table([a; b_coeff'], 'VariableNames', {'Value'}, ...
'RowNames', {'a', 'b', 'c', 'd', 'e', 'f'});
r2_table = table(R2, 'VariableNames', {'R2_Value'});
filename = 'matlab_nonlinear_regression1.xlsx';
writetable(results_table, filename, 'Sheet', 'Storage Comparison');
writetable(coeff_table, filename, 'Sheet', 'Regression Coefficients', 'WriteRowNames', true);
writetable(r2_table, filename, 'Sheet', 'R2 Value');
disp('Excel report generated successfully!');
figure;
scatter(S, S_Estimate, 'b', 'filled');
hold on;
plot([min(S), max(S)], [min(S), max(S)], 'r--', 'LineWidth', 2);
hold off;
xlabel('Actual Storage (S)');
ylabel('Estimated Storage (S_{Estimate})');
title(['Nonlinear Regression (R^2 = ', num2str(R2), ')']);
legend('Data Points', '1:1 Line', 'Location', 'Best');
grid on;

类别

在 帮助中心 和 File Exchange 中查找有关 Linear and Nonlinear Regression 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by