fmincon nonlinear inequality constraints

4 次查看(过去 30 天)
Hello,
I am working on developing deterioration model for road pavement condition data for my thesis. I have been using the Matlab curve fitting tool with a custom deterioration model. I have my x (age of pavement) and y (condition index) data. my model equation is y = o-exp(a-b*c^log(1/x)). I am estimating coefficient a, b, c and o based on some constraints. I was able to do all that but i am stuck with a constrain that I dont know how to formulate.
I just learn about fmincon, so i decided to use fmincon to optimized my parameter a, b, c and o while minimizing the SSE.
My main issue is how to set up the nonlinear inequality constraint so that:
function [c, ceq]:
if x = 2.7
c = o-exp(a-b*c^log(1/x)) is <= 35
ceq =[ ]
end
My main issue is that when i use a, b, c and o in my inequality c, the model does not work. but when i use number value for a, b, c and o my model work. but it seems like it doesnt follow my constraint.
I will appreciate if anyone can help me with this issue. Thank you in advance.
below is my code and attached the excel data. Thank you in advance and let me know if you have any question.
function [yp,objective,c,popt] = regX(xm, ym)
% clear session, close plots, clear screen
clear all; close all; clc
% data for regression
S = uiimport('-file'); %this gives option to automatically select excel file
% It will show an error message but wait until a pop window show
xm = S.data(:,1);
ym = S.data(:,2);
[xData, yData] = prepareCurveData( xm, ym );
% define prediction function
%yp = @(p) p(1) + p(2)./xm + p(3).*log(xm);
yp = @(p) p(4)-exp(p(1)-p(2).*p(3).^log(1./xm))
% initial parameter guess
p0 = [2.964, 3.825, 1.231, 5]; %pdi con GR
%p0 = [25.706, 26.482, 1.02, 4.6]; %psi con & crc GR
% define objective function (scaled sum of squared errors)
%objective = @(p) sum(((yp(p)-ym)./ym).^2);
objective = @(p) sum(((yp(p)-ym)./ym).^2);
disp(['Initial Objective: ' num2str(objective(p0))])
% linear constraints
A = [];
b = [];
Aeq = [];
beq = [];
% variable bounds
lb = []; % ones(3)*0.2;
ub = []; % ones(3)*1.5;
% create file nlcon.m for nonlinear constraints
function [c,ceq] = nlcon(xm)
for xm = 10
c = (5-exp(2.964-3.825.*1.231.^log(1./xm))) - 2.7; %pdi CON GR
%c = (4.6-exp(25.706-26.482.*1.02.^log(1./xm)) - 2.7); %psi con & crc GR
%c = (p(4)-exp(p(1)-p(2).*p(3).^log(1./xm)) - 2.7); % This is the
%inequality i am trying to use with parameter to be estimated in
%the regression
%ceq = sum(x.^2) - 40;
ceq = [];
end
end
%nonlincon = @nlcon;
% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]
% = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
%options = optimoptions('fmincon', 'Algorithm', 'sqp', 'Display', 'iter-detailed',...
%'MaxFunctionEvaluations', 100000, 'MaxIterations', 2000,...
%'FunctionTolerance', 1e-10);
popt = fmincon(objective,p0,A,b,Aeq,beq,lb,ub,@(xm) nlcon(xm));
% print results
disp(['Final Objective: ' num2str(objective(popt))])
disp(['Optimal parameters: ' num2str(popt)])
% plot results
plot(xm,ym,'ro')
hold on
plot(xm,yp(p0),'bx')
plot(xm,yp(popt),'gs')
legend('measured','initial predicted','optimal predicted')
ylabel('y')
xlabel('x')
end

采纳的回答

Matt J
Matt J 2020-11-15
function [c, ceq] = nlcon(p)
c = p(4)-exp(p(1)-p(2).*p(3).^log(1./2.7)) -35 ;
ceq =[ ];
end
  6 个评论
Amara Kouyate
Amara Kouyate 2020-11-16
okay thanks Matt. I thought that i was doing something worng. Do you know any function or anything that i can do so that my solution won't stuck at local minimum? Thank you!
Amara Kouyate
Amara Kouyate 2020-11-17
hello Matt,
I was able to fix the issue of having same value with fmincon. I was using the worng input with the @nlcon function. Thank you very much for your help! if you get a chance, don't about my previous question. This the correct formulation of the @nlcon:
function [c, ceq] = nlcon(p)
c = p(4)-exp(p(1)-p(2).*p(3).^log(1./35)) -2.7 ;
ceq =[ ];
end

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Systems of Nonlinear Equations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by