Optimisation problem in matlab

Hello guys,
So I have this scatte plot in the form of a measured impedance according to frequency Z(f), and I have an RL ladder circuit seen in the following figure
.
I need some help in writing an optimisation problem that will allow me to get the values of the parameters R_i L_i that will eventually give an equivalent impedance of the RL ladder similar to the one measured- Z(f) -.
Thanks in advance.
Wallflower

2 件のコメント

Rik
Rik 2021 年 6 月 21 日
This looks like a homework assignment. You can find guidelines for posting homework on this forum here. If you have trouble with Matlab basics you may consider doing the Onramp tutorial (which is provided for free by Mathworks). If your main issue is with understanding the underlying concept, you may consider re-reading the material you teacher provided and ask them for further clarification.
wallflower
wallflower 2021 年 6 月 21 日
It is not a homework assignment. I am trying to reproduce a methology I've seen in a scientific paper.

サインインしてコメントする。

回答 (1 件)

Sergey Kasyanov
Sergey Kasyanov 2021 年 6 月 21 日
編集済み: Torsten 2022 年 7 月 12 日

0 投票

Hello!
Try that:
% tune only next 2 lines
L = 5;% steps in ladder
Z_max = 1e2;% max Z for one R or Xl
% goal Z(w)
% w0 - array with measured angular frequences
% Z0 - array with measured impedances
% R0 - Rs_1.txt
% L0 - 2*Self_Energy_1.txt
w0 = (2.*pi.*[[1e-3,50,500],1e3:2e3:500e3])';
% Start inserted to make code work
R0 = rand(size(w0));
L0 = rand(size(w0));
% End inserted to make code work
Z0 = R0(:,1) + 1i*w0.*L0(:,1);
w0 = reshape(w0, 1, []);
Z0 = reshape(Z0, 1, []);
%% calculate z(w) for ladder
R = sym('R', [1,L+1]);
L = sym('L', [1,length(R)]);
w = sym('w');
par = @(a, b) 1/(1/a + 1/b);
h1 = @(a, i) R(i) + par(a, 1i*w*L(i));
Z = R(end);
for i = (length(R)-1):-1:1
Z = h1(Z, i);
end
% get lambda function for optimization
fun1 = matlabFunction(Z);
% do not look at that code
s = 'fun2 = @(vals, w) fun1(';
for i = 1:(length(R)*2 - 1)
s = [s, sprintf('vals(%i), ', i)];
end
s = [s, 'w);'];
eval(s);
% create optimization function
fun3 = @(vals) sum(abs(fun2(vals, w0) - Z0).^2);
% number of variables
N = length(symvar(Z)) - 1;
%the better way is to use more powerfull algorithms such as genetic algoritms which can stands
%against a huge amount of local min. There are a lot of algorithms that
%realized in matlab toolboxes. Try!
options = optimoptions('ga', 'Display', 'iter', 'PopulationSize', N*200, 'MaxGenerations', N*200);
res = ga(fun3, N, [], [], [], [], zeros(1, N), Z_max*[ones(1,floor(N/2))/2/pi, ones(1, ceil(N/2))], [], options);
Single objective optimization: 11 Variable(s) Options: CreationFcn: @gacreationuniform CrossoverFcn: @crossoverscattered SelectionFcn: @selectionstochunif MutationFcn: @mutationadaptfeasible Best Mean Stall Generation Func-count f(x) f(x) Generations 1 4400 2.953e+14 2.953e+14 0 2 6490 2.953e+14 2.953e+14 0 3 8580 2.953e+14 2.953e+14 0 4 10670 2.953e+14 2.953e+14 0 5 12760 2.953e+14 2.953e+14 0 6 14850 2.953e+14 2.953e+14 0 7 16940 2.953e+14 2.953e+14 0 8 19030 2.953e+14 2.953e+14 0 9 21120 2.953e+14 2.953e+14 0 10 23210 2.953e+14 2.953e+14 0 11 25300 2.953e+14 2.953e+14 0 12 27390 2.953e+14 2.953e+14 0 13 29480 2.953e+14 2.953e+14 0 14 31570 2.953e+14 2.953e+14 0 15 33660 2.953e+14 2.953e+14 0 16 35750 2.953e+14 2.953e+14 0 17 37840 2.953e+14 2.953e+14 0 18 39930 2.953e+14 2.953e+14 0 19 42020 2.953e+14 2.953e+14 0 20 44110 2.953e+14 2.953e+14 0 21 46200 2.953e+14 2.953e+14 0 22 48290 2.953e+14 2.953e+14 0 23 50380 2.953e+14 2.953e+14 0 24 52470 2.953e+14 2.953e+14 1 25 54560 2.953e+14 2.953e+14 0 26 56650 2.953e+14 2.953e+14 0 27 58740 2.953e+14 2.953e+14 0 28 60830 2.953e+14 2.953e+14 0 29 62920 2.953e+14 2.953e+14 0 30 65010 2.953e+14 2.953e+14 1 Best Mean Stall Generation Func-count f(x) f(x) Generations 31 67100 2.953e+14 2.953e+14 0 32 69190 2.953e+14 2.953e+14 1 33 71280 2.953e+14 2.953e+14 2 34 73370 2.953e+14 2.953e+14 0 35 75460 2.953e+14 2.953e+14 0 36 77550 2.953e+14 2.953e+14 0 37 79640 2.953e+14 2.953e+14 0 38 81730 2.953e+14 2.953e+14 0 39 83820 2.953e+14 2.953e+14 1 40 85910 2.953e+14 2.953e+14 2 41 88000 2.953e+14 2.953e+14 3 42 90090 2.953e+14 2.953e+14 0 43 92180 2.953e+14 2.953e+14 0 44 94270 2.953e+14 2.953e+14 0 45 96360 2.953e+14 2.953e+14 0 46 98450 2.953e+14 2.953e+14 1 47 100540 2.953e+14 2.953e+14 2 48 102630 2.953e+14 2.953e+14 3 49 104720 2.953e+14 2.953e+14 4 50 106810 2.953e+14 2.953e+14 0 51 108900 2.953e+14 2.953e+14 0 52 110990 2.953e+14 2.953e+14 0 53 113080 2.953e+14 2.953e+14 1 54 115170 2.953e+14 2.953e+14 0 55 117260 2.953e+14 2.953e+14 1 Optimization terminated: average change in the fitness value less than options.FunctionTolerance.
% first floor(N/2) of res - L(1), L(2), ...
% last ceil(N/2) of res - R(1), R(2), ...
%plot results
figure;plot(w0, abs(Z0), w0, abs(fun2(res,w0)));
figure;plot(w0, angle(Z0), w0, angle(fun2(res,w0)));

14 件のコメント

wallflower
wallflower 2021 年 6 月 21 日
Thanks for your response, but I get the following error:
Error using fminbnd (line 109)
Bounds must be scalars.
Sergey Kasyanov
Sergey Kasyanov 2021 年 6 月 21 日
sorry, use fmincon instead:
res = fmincon(fun,rand(1,N),[],[],[],[],zeros(1,N),9e9*ones(1,N),[],options);
wallflower
wallflower 2021 年 6 月 21 日
No worries.
I tried it and the fit and the measurement are no where near each other x(, the simulation only found a local minimum ...
Sergey Kasyanov
Sergey Kasyanov 2021 年 6 月 22 日
What is about ga()?
wallflower
wallflower 2021 年 6 月 22 日
For the ga() this is the result I get :'(
Sergey Kasyanov
Sergey Kasyanov 2021 年 6 月 22 日
Can you load Z(f) here as .mat file?
wallflower
wallflower 2021 年 6 月 22 日
Yes sure. Please find attachated the Z(f) file (a frequency-dependent resistance whose values have been computed using finite element method to take into account skin & proximity effects on the following angular frequency range of interest: w0 = 2.*pi.*[[1e-3,50,500],1e3:2e3:500e3]
Sergey Kasyanov
Sergey Kasyanov 2021 年 6 月 24 日
Are Z(f) pure resistance? That means that it does not depend on frequency, so problem is not correct.
I corrected code in answer.
wallflower
wallflower 2021 年 6 月 24 日
Z(f) is pure resistance and it does depend on frequency because of skin and proximity effects. But I've had time to search a little and I can do the same ladder circuit on a complex impedance by including the inductive part of the cable so now Z(f) becomes: R(f)+j.*L(f).*w0. I'll try your code on this complex impedance. Attached are the files that I will be using.
Note: L(f) is deduced from the self_energy file as follows: L(f) = 2*Self_energy_values
Sergey Kasyanov
Sergey Kasyanov 2021 年 6 月 24 日
I corrected code. There were some misprints.
Now it is working, but It is very sensitive to limits of Z_max and amount of ladder steps. Try, please. Does it working on your machine?
Anes MESSADI
Anes MESSADI 2022 年 7 月 12 日
I tried to use the code you proposed but it doesn't work, can you please share the corrected one.
Thank you,
Sergey Kasyanov
Sergey Kasyanov 2022 年 7 月 12 日
I can't test and correct code at the moment. What error is there?
Anes MESSADI
Anes MESSADI 2022 年 7 月 12 日
there are no errors, it works, but the fitting is not good, I used only the Rs_1.txt and Self_Energy_1.txt just to test the code.
Thank you for your help
Sergey Kasyanov
Sergey Kasyanov 2022 年 7 月 13 日
Try to use that code with your w0 and Z0. Also you need function ladder_z in the file.
% tune only next 2 lines
L = 10;% steps in ladder
Z_max = 1e2;% max Z for one R or Xl
% w0 - array with measured angular frequences
% Z0 - array with measured impedances
% create data for testing
w0 = 2 * pi * logspace(-3, 2, 10);
z0 = [rand(1, L);rand(1, L)];
Z0 = [];
for i = 1:length(w0)
Z0(i) = ladder_z(z0(1, :) + z0(2, :) * 1i * w0(i));
end
% main code
h1 = @(vals, w) ladder_z(vals(1:fix(length(vals) / 2)) + 1i * w * vals(fix(length(vals) / 2) + 1:end));
h2 = @(vals) arrayfun(@(iw) h1(vals, iw), w0);
f = @(vals) sum(abs(h2(vals) - Z0).^2);
options = optimoptions('ga', 'Display', 'iter', 'PopulationSize', L*200, 'MaxGenerations', L*200, 'UseParallel', true);
res = ga(f, L * 2, [], [], [], [], zeros(1, 2 * L), Z_max*[ones(1, L), ones(1, L) / 2 / pi], [], options);
% check
figure;plot(w0, abs(Z0), w0, abs(h2(res)));

サインインしてコメントする。

カテゴリ

ヘルプ センター および File ExchangeLadder Diagram Integration についてさらに検索

製品

リリース

R2019b

質問済み:

2021 年 6 月 21 日

コメント済み:

2022 年 7 月 13 日

Community Treasure Hunt

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

Start Hunting!

Translated by