Could anybody help me with an iterative process?

Hello!

First, I am sorry if my english is not good. I will try to make the doubt as clear as possible.

The code shown below is a function that is suppose to iterate some areas (stored in B0). B0 is a matrix 27x6 in this case, where each row represents a different system with 6 different areas. What I want to do is to minimize my values of B0 (maintaining their minimum value to 1e-4), while making RF equal to 1, with a tolerance of 0.05 (i.e. 0.95<= RF <= 1.05 would be OK). The thing is, I am not able to make it converge, and moreover, sometimes I get negative values for B0, which would not make sense.

I think I made the code a little messy, so I don't know if I am explaining myself properly, but if anybody is willing to help me, I would really appreciate it. Thank you!!

P.D: y,x,Mx are fixed. Y and X are vector positions (27x6, 1 row for each system, 6 columns, each one describing the position of each area), while Mx is the value of a moment.

idx = size(y);
idx = idx(2);
B0=ones(length(y),idx);
B0= B0.*0.001; % Initial guess for the area of the booms ~ 100 mm^2
correct=0;
Bmin=1e-4;
sigma_max = sigma_max*0.7;
    a = sum(B0,2);
    CGx = sum(((B0.*x)./a),2);
    I = sum(B0.*y.^2,2);
    sigma = (y.*Mx)./I;
    RF = abs(sigma_max./sigma);
    k1 = find(RF<0.95);
    k2 = find(RF>1.05);
    k3 = find(1.05>RF & RF>0.95);
    k4 = find(B0<Bmin);
    count = 0;
    mySize = size(B0);
    tolRF = 0.05;
    while (~isempty(k1) || ~isempty(k2) || ~isempty(k4))
        count = count+1;
        if(~isempty(k1))
        B0(k1) = B0(k1)+0.00001;
        end
        if(~isempty(k2))
        B0(k2) = B0(k2)-0.00001;
        end
        if(~isempty(k3))
        B0(k3) = B0(k3);
        end
        a = find(B0<Bmin);
        if(~isempty(a))
            B0(a) = Bmin;
        end
        k1 = find(RF<(1-tolRF));
        k2 = find(RF>(1+tolRF));
        k3 = find((1+tolRF)>RF & RF>(1-tolRF));
        k4 = a;
        I = sum(B0.*y.^2,2);
        sigma = abs((y.*Mx)./I);
        RF = (sigma_max./sigma);
    end
end

Answer 1

[Chetan]

Hello Juan,

I see that you intend to refine the `B0` matrix so that the stress ratio `RF` is constrained within 0.95 to 1.05, while also ensuring that the values of `B0` do not fall below the threshold of `1e-4`.

There are issues in the current code that might prevent it from converging and could lead to `B0` taking negative values.

Here are the steps to rectify the problem:

• Guarantee that `B0` doesn't decrease beneath `Bmin`.
• Implement proportional changes to `B0` contingent on the distance of `RF` from the target range.
• Re-evaluate `RF` after every modification to `B0`.
• Ensure `sigma_max` is accurately defined and incorporated.

Below is the modified function:

% Example inputs for y, x, Mx, and sigma_max
y = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6]; % Sample y coordinates
x = [0.5, 0.4, 0.3, 0.2, 0.1, 0];   % Sample x coordinates
Mx = 100;                            % Sample bending moment
sigma_max = 200;                     % Sample maximum stress

% Convert y and x to column vectors if they aren't already
y = y(:);
x = x(:);

% Execute the function with the sample inputs
optimized_B0 = optimize_areas(y, x, Mx, sigma_max);

Convergence achieved.
Iterations: 887

% Output the optimized boom areas
disp('Optimized B0:');

Optimized B0:

disp(optimized_B0);

    6.8086
    3.4268
    2.2788
    1.7076
    1.3583
    1.1356

% The function to optimize boom areas
function B0 = optimize_areas(y, x, Mx, sigma_max)
% Define minimum boom area, tolerance, and initial guess
Bmin = 1e-4;
tolRF = 0.05;
B0 = ones(size(y)) * 0.001;
sigma_max = sigma_max * 0.7; % Modify sigma_max as necessary

% Initial calculations for area, centroid, and stress
a = sum(B0, 2);
I = sum(B0 .* y.^2, 2);
sigma = (y .* Mx) ./ I;
RF = abs(sigma_max ./ sigma);

% Iteration counter and limit
count = 0;
maxIterations = 10000;

% Iterative process to adjust B0
while any(RF < (1 - tolRF) | RF > (1 + tolRF)) && count < maxIterations
    count = count + 1;
    for i = 1:numel(B0)
        if RF(i) < (1 - tolRF)
            B0(i) = B0(i) + B0(i) * 0.01; % Increment by 1%
        elseif RF(i) > (1 + tolRF)
            B0(i) = B0(i) - B0(i) * 0.01; % Decrement by 1%
        end
        B0(i) = max(B0(i), Bmin); % Enforce minimum value
    end
    
    % Recompute values after updating B0
    a = sum(B0, 2);
    I = sum(B0 .* y.^2, 2);
    sigma = (y .* Mx) ./ I;
    RF = abs(sigma_max ./ sigma);
end

% Check for convergence
if count == maxIterations
    disp('Maximum iterations reached without convergence.');
else
    disp('Convergence achieved.');
end
disp(['Iterations: ', num2str(count)]);
end

Hope this should help resolve the issues

Best regards,

Chetan