different results with fmincon

I would like to ask about the relationship between the objective function and nonlinear constraint function defined in the fmincon function, as they produce different results. when I write code as follows

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[10,10,10,10,10,10];
B=0.8;
v=30/3.6;
x = fmincon(@(x)fun(x,B,v),x0,[],[],[],[],lb,ub,@(x)nonlcon(x,B,v))
a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);
function [f,difference]=fun(x,B,v)% 'difference' is defined
s = tf('s');
w_values = 0:0.5:100;
G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
[mag] = bode(G, w_values);
k_values = zeros(size(w_values));
difference=zeros(size(w_values));%  The 'difference' is initialized to array 0
for i = 1:length(w_values)
    k = exp(-w_values(i)*B/v);
    k_values(i) = k; 
end
for i = 1:length(w_values)
   difference(i)=abs(mag(i)-k_values(i));% calculate 'difference'
end
f=sum(difference);
end
function [c,ceq] = nonlcon(x,B,v)
[difference]=fun(x,B,v);
%Is the 'difference' here referring to the difference in the function [f,difference]=fun(x,B,v)?"
c=max(difference)-0.001;
ceq=[];
end

results:x =

9.9883 0.0000 0.0021 10.0000 2.6715 0.1221

Directly using "difference" from the objective function as an input parameter for the nonlinear constraint function resulted in different outcomes. Is it reasonable to reference "difference" directly in this way? Do I need to call it like [difference] = fun(x, B, v); as shown in the previous code snippet?

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[10,10,10,10,10,10];
x = fmincon(@fun,x0,[],[],[],[],lb,ub,@nonlcon)%Here it has been altered
a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);
function [f]=fun(x)% Here it has been altered
s = tf('s');
B=0.8;
v=30/3.6;
w_values = 0:0.5:100;
G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
[mag] = bode(G, w_values);
k_values = zeros(size(w_values));
difference=zeros(size(w_values));
for i = 1:length(w_values)
    k = exp(-w_values(i)*B/v);
    k_values(i) = k; 
end
for i = 1:length(w_values)
   difference(i)=abs(mag(i)-k_values(i));
end
f=sum(difference);
end
function [c,ceq] = nonlcon(difference)%Here it has been altered
c=max(difference)-0.001;
ceq=[];
end

results=

1.0e-03 *

0.7235 0.0006 0.0001 0.8284 0.1887 0.0093

Answer 1

[Matt J]

function [c,ceq] = nonlcon(x,B,v)

  [~,difference]=fun(x,B,v);

  c=max(difference)-0.001;
  ceq=[];
end

[guo qing]

Thanks you very much.Can the 'difference' variable be directly defined as an input parameter for the nonlinear constraint? Will the 'difference' be passed from the objective function to the nonlinear constraint function in this way?

[guo qing]

And I just tried the modification method you mentioned, [~,difference]=fun(x,B,v); The final fitting effect is not as good as before

[Walter Roberson]

The above would have to be combined with

function [f, difference]=fun(x)% Here it has been altered

[guo qing]

Sorry, I didn't quite understand. Could you provide the complete code for that?

[Walter Roberson]

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[10,10,10,10,10,10];
x = fmincon(@fun,x0,[],[],[],[],lb,ub,@nonlcon)%Here it has been altered
a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);

function [f, difference]=fun(x)% Here it has been altered
    s = tf('s');
    B=0.8;
    v=30/3.6;
    w_values = 0:0.5:100;
    G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
    [mag] = bode(G, w_values);
    k_values = zeros(size(w_values));
    difference=zeros(size(w_values));
    for i = 1:length(w_values)
        k = exp(-w_values(i)*B/v);
        k_values(i) = k; 
    end
    for i = 1:length(w_values)
       difference(i)=abs(mag(i)-k_values(i));
    end
    f=sum(difference);
end

function [c,ceq] = nonlcon(difference)%Here it has been altered
    [~,difference] = fun(x,B,v);
    c = max(difference)-0.001;
    ceq = [];
end

[guo qing]

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[10,10,10,10,10,10];
x = fmincon(@fun,x0,[],[],[],[],lb,ub,@nonlcon)%Here it has been altered

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

x = 1×6

1.0e-03 * 0.6350 0.0004 0.0001 0.7264 0.1657 0.0082

a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);

function [f]=fun(x)% Here it has been altered
    s = tf('s');
    B=0.8;
    v=30/3.6;
    w_values = 0:0.5:100;
    G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
    [mag] = bode(G, w_values);
    k_values = zeros(size(w_values));
    difference=zeros(size(w_values));
    for i = 1:length(w_values)
        k = exp(-w_values(i)*B/v);
        k_values(i) = k; 
    end
    for i = 1:length(w_values)
       difference(i)=abs(mag(i)-k_values(i));
    end
    f=sum(difference);
end

function [c,ceq] = nonlcon(difference)%Here it has been altered
    c = max(difference)-0.001;
    ceq = [];
end

However, it works to directly use 'difference' as an input variable without calling the 'fun' function, and the final results are also very good. Why is this the case

[Matt J]

However, it works to directly use 'difference' as an input variable without calling the 'fun' function

fmincon will always pass the unknown parameter vector x to nonlcon, regardless of whether you name the input difference or something else. Therefore what you have implemented ciould easily have been rewritten as

function [c,ceq] = nonlcon(x)
    c = max(x)-0.001;
    ceq = [];
end

Moreover, this constraint is not nonlinear at all and is equivalent to simply tightening the upper bound to ub=0.001*ones(1,6).

As for why you see results you like better, it would seem that the solution you need really lies within the narrowwe region 0<=x(i)<=0.001. By tightening ub to this region, you have made that solution easier to find.

[Torsten]

I think what you really wanted to implement was

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[10,10,10,10,10,10];
x = fmincon(@fun,x0,[],[],[],[],lb,ub,@nonlcon,optimoptions('fmincon','Algorithm','interior-point'))%Here it has been altered

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.


Consider enabling the interior point method feasibility mode.

x = 1×6

3.0356 8.6287 0.0157 3.2221 9.4031 2.4835

a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);

function [f,difference]=fun(x)% Here it has been altered
    s = tf('s');
    B=0.8;
    v=30/3.6;
    w_values = 0:0.5:100;
    G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
    [mag] = bode(G, w_values);
    k_values = zeros(size(w_values));
    difference=zeros(size(w_values));
    for i = 1:length(w_values)
        k = exp(-w_values(i)*B/v);
        k_values(i) = k; 
    end
    for i = 1:length(w_values)
       difference(i)=abs(mag(i)-k_values(i));
    end
    f=sum(difference);
end

function [c,ceq] = nonlcon(x)%Here it has been altered
    [~,difference] = fun(x);
    c = max(difference)-0.001;
    ceq = [];
end

Note that G only has 5 independent parameters instead of 6 (you can always normalize x(3) or x(6) to 1). This makes your problem singular.

[Sam Chak]

Hey @guo qing, could you figure out which system produces the desired frequency response? or ?

s   = tf('s');
x0  = [1,1,1,1,1,1];
lb  = [0,0,0,0,0,0];
ub  = [10,10,10,10,10,10];

x1  = fmincon(@fun,x0,[],[],[],[],lb,ub,@nonlcon)%Here it has been altered

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

x1 = 1x6

1.0e-03 * 0.6350 0.0004 0.0001 0.7264 0.1657 0.0082

G1  = (x1(1) + x1(2)*s + x1(3)*s^2) / (x1(4) + x1(5)*s + x1(6)*s^2)

G1 =
 
  1.265e-07 s^2 + 4.499e-07 s + 0.000635
  ---------------------------------------
  8.168e-06 s^2 + 0.0001657 s + 0.0007264
 
Continuous-time transfer function.

bode(G1), grid on
title('Bode Diagram of System G_1')

x2  = fmincon(@fun2,x0,[],[],[],[],lb,ub,@nonlcon2,optimoptions('fmincon','Algorithm','interior-point'))%Here it has been altered

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.


Consider enabling the interior point method feasibility mode.

x2 = 1x6

3.0356 8.6287 0.0157 3.2221 9.4031 2.4835

G2  = (x2(1) + x2(2)*s + x2(3)*s^2) / (x2(4) + x2(5)*s + x2(6)*s^2)

G2 =
 
  0.01569 s^2 + 8.629 s + 3.036
  -----------------------------
   2.484 s^2 + 9.403 s + 3.222
 
Continuous-time transfer function.

bode(G2), grid on
title('Bode Diagram of System G_2')

function [f]=fun(x)% Here it has been altered
    s = tf('s');
    B=0.8;
    v=30/3.6;
    w_values = 0:0.5:100;
    G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
    [mag] = bode(G, w_values);
    k_values = zeros(size(w_values));
    difference=zeros(size(w_values));
    for i = 1:length(w_values)
        k = exp(-w_values(i)*B/v);
        k_values(i) = k; 
    end
    for i = 1:length(w_values)
       difference(i)=abs(mag(i)-k_values(i));
    end
    f=sum(difference);
end

function [c,ceq] = nonlcon(difference)%Here it has been altered
    c = max(difference)-0.001;
    ceq = [];
end

function [f,difference]=fun2(x)% Here it has been altered
    s = tf('s');
    B=0.8;
    v=30/3.6;
    w_values = 0:0.5:100;
    G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
    [mag] = bode(G, w_values);
    k_values = zeros(size(w_values));
    difference=zeros(size(w_values));
    for i = 1:length(w_values)
        k = exp(-w_values(i)*B/v);
        k_values(i) = k; 
    end
    for i = 1:length(w_values)
       difference(i)=abs(mag(i)-k_values(i));
    end
    f=sum(difference);
end

function [c,ceq] = nonlcon2(x)%Here it has been altered
    [~,difference] = fun2(x);
    c = max(difference)-0.001;
    ceq = [];
end

[guo qing]

@Matt J Thank you very much.Now I know that difference is actually the original x defined, but the interesting thing is that I don't pass the 'difference' defined by the 'fun' function to 'nonclon', as you said just set the upper bound of x to 0.001 to get a better fit. Instead, I passed 'difference' to 'nonclon' and the fitting effect got worse. This is my source file and drawing file. Do you have any suggestions on how to get better fitting effect

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[0.001,0.001,0.001,0.001,0.001,0.001];
B=0.8;
v=30/3.6;
options = optimoptions('fmincon', 'MaxIterations', 5000);
x = fmincon(@(x)fun(x,B,v),x0,[],[],[],[],lb,ub,[], options)
a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);
function [f,difference]=fun(x,B,v)
s = tf('s');
w_values = 0:0.5:100;
G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
[mag,~]= bode(G, w_values);
k_values = zeros(size(w_values));
difference=zeros(size(w_values));% 初始化存储 k 值的数组
for i = 1:length(w_values)
    k = exp(-w_values(i)*B/v);
    k_values(i) = k; % 将每次计算得到的 k 值存储在数组中
end
for i = 1:length(w_values)
   difference(i)=abs(mag(i)-k_values(i));
end
f=sum(difference);
end

plot file as following

B = 0.8;
v = 30/3.6;
w_values = 0:0.5:100;
new_function = exp(-w_values*B/v);
% 计算拟合传递函数的幅度响应
x = [x(1), x(2), x(3), x(4), x(5), x(6)];
s = tf('s');
G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
[mag,~]= bode(G, w_values);
mag = [squeeze(mag)]';
% 计算差值
k = abs(new_function - mag);
max(k)
% 绘制拟合传递函数和给定函数的幅频特性曲线
figure;
plot(w_values,mag,'LineWidth',1.5)
hold on
plot(w_values,new_function,'LineWidth', 1.5);
hold on
title('幅频特性比较');
xlabel('角频率(rad/s)');
ylabel('幅值');
legend('拟合模型', '参照模型');
grid on;
ylim([-0.1,1.01]);%设置纵坐标范围
xticks = 0:20:100;
set(gca, 'XTick', xticks)
% 绘制差值与频率的图像
figure;
plot(w_values,k,'LineWidth',1);
hold on
plot(w_values,zeros(size(w_values)),'r--');
xlabel('角频率(rad/s)');
ylabel('拟合误差');
grid on;
ylim([-0.025,max(k)+0.01]); % 设置纵坐标范围
xticks = 0:20:100;
set(gca, 'XTick', xticks)

[guo qing]

@Sam Chak Thank you very much.I want to get a good function fitting effect.This is my source file and drawing file

x0=[1,1,1,1,1,1];
lb=[0,0,0,0,0,0];
ub=[0.001,0.001,0.001,0.001,0.001,0.001];
B=0.8;
v=30/3.6;
options = optimoptions('fmincon', 'MaxIterations', 5000);
x = fmincon(@(x)fun(x,B,v),x0,[],[],[],[],lb,ub,[], options)
a0 = x(1);
a1 = x(2);
a2 = x(3);
b0 = x(4);
b1 = x(5);
b2 = x(6);
function [f,difference]=fun(x,B,v)
s = tf('s');
w_values = 0:0.5:100;
G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
[mag,~]= bode(G, w_values);
k_values = zeros(size(w_values));
difference=zeros(size(w_values));% 初始化存储 k 值的数组
for i = 1:length(w_values)
    k = exp(-w_values(i)*B/v);
    k_values(i) = k; % 将每次计算得到的 k 值存储在数组中
end
for i = 1:length(w_values)
   difference(i)=abs(mag(i)-k_values(i));
end
f=sum(difference);
end

plot file

B = 0.8;
v = 30/3.6;
w_values = 0:0.5:100;
new_function = exp(-w_values*B/v);
% 计算拟合传递函数的幅度响应
x = [x(1), x(2), x(3), x(4), x(5), x(6)];
s = tf('s');
G = (x(1) + x(2)*s + x(3)*s^2) / (x(4) + x(5)*s + x(6)*s^2);
[mag,~]= bode(G, w_values);
mag = [squeeze(mag)]';
% 计算差值
k = abs(new_function - mag);
max(k)
% 绘制拟合传递函数和给定函数的幅频特性曲线
figure;
plot(w_values,mag,'LineWidth',1.5)
hold on
plot(w_values,new_function,'LineWidth', 1.5);
hold on
title('幅频特性比较');
xlabel('角频率(rad/s)');
ylabel('幅值');
legend('拟合模型', '参照模型');
grid on;
ylim([-0.1,1.01]);%设置纵坐标范围
xticks = 0:20:100;
set(gca, 'XTick', xticks)
% 绘制差值与频率的图像
figure;
plot(w_values,k,'LineWidth',1);
hold on
plot(w_values,zeros(size(w_values)),'r--');
xlabel('角频率(rad/s)');
ylabel('拟合误差');
grid on;
ylim([-0.025,max(k)+0.01]); % 设置纵坐标范围
xticks = 0:20:100;
set(gca, 'XTick', xticks)

[Matt J]

Do you have any suggestions on how to get better fitting effect

We already know how to get a better fit. Set ub=0.001.

[guo qing]

Thank you very much. It means that optimizing from the perspective of the total sum of differences is better than setting each individual component of the difference to be less than a certain value. Best wishes.